The spin expansion for binary black hole mergers: 
new predictions and future directions 
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In a recent paper we introduced a spin expansion that provides a simple yet powerful way 
to understand aspects of binary black hole (BBH) merger. This approach relies on the symmetry 
properties of initial and final quantities like the black hole mass m, kick velocity k, and spin vector 
s, rather than a detailed understanding of the merger dynamics. In this paper, we expand on this 
proposal, examine how well its predictions agree with current simulations, and discuss several future 
directions that would make it an even more valuable tool. The spin expansion yields many new 
predictions, including several exact results that may be useful for testing numerical codes. Some of 
these predictions have already been confirmed, while others await future simulations. We explain 
how a relatively small number of simulations — 10 equal-mass simulations, and 16 unequal-mass 
simulations — may be used to calibrate all of the coefficients in the spin expansion up to second 
order at the minimum computational cost. For a more general set of simulations of given covariance, 
we derive the minimum-variance unbiased estimators for the spin expansion coefficients. We discuss 
how this calibration would be interesting and fruitful for general relativity and astrophysics. Finally, 
we sketch the extension to eccentric orbits. 
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I. INTRODUCTION 

Binary black hole (BBH) merger — in which two spin- 
ning black holes inspiral due to the emission of gravi- 
tational radiation and eventually merge to form a sin- 
gle spinning black hole — is one of the most impor- 
tant problems in classical general relativity, and has sig- 
nificant ramifications for astrophysics, cosmology, and 
gravitational-wave observations. For decades, analytical 
and numerical approaches to the BBH merger problem 
have been frustrated by conceptual and technical difficul- 
ties associated with solving the non-linear Einstein equa- 
tions — especially during the last few orbits and final 
plunge, when the "luminosity" in gravitational radiation 
is highest. Dramatic progress came in 2005, as new in- 
sights and increased computational resources finally al- 
lowed numerical relativists to simulate the entire merger 
— including the last few orbits of inspiral, the plunge, 
the formation of a common event horizon, and the ring- 
down of the final Kerr black hole 0, d, Q ■ Following this 
breakthrough, simulations of BBH merger have produced 
a number of remarkable results. 

The most surprising and interesting results have been 
obtained just within the past year, from simulations of 
mergin g b lack holes with larg e initial spins [a, lO , III, [a , l 9| , 

[13, [ni, III, m, H, M, [n, eje m, Hfl, minii [11,121 . 

The result that has received the most attention is that 
highly-spinning initial black holes can merge to form a 
final black hole with an enormous recoil velocity — as 
large as 4, 000 km/s — relative to the binary's center-of- 
momentum frame [E, [3, [3 ■ The idea of a supermassive 
black hole rocketing through its host galaxy at such a 
speed has understandably caught the attention of many 
astrophysicists! 

In this paper, we highlight a second surprising result. 
Despite the complicated and non-linear dynamics of the 
merger process, the final state of the merger seems, in 



some sense, to be an unexpectedly simple and smooth 
function of the initial state of the binary. We would like 
to make this statement more quantitative and precise. 

In a recent paper [l| , we introduced a "spin expansion" 
formalism for understanding aspects of BBH merger. We 
expand on this proposal in several ways in the present 
paper. Here is a recap of the basic idea. Even though 
the merger is a messy non-linear process, it is useful to 
regard it as a map from a simple initial state (two well 
separated Kerr black holes with mass ratio q = Mb/ Ma 
and dimensionless spins a and b) to a simple final state 
(a final Kerr black with mass m, spin vector s and kick 
velocity k). Given any final quantity / (e.g. m, k, or 
s), we can Taylor expand the function /(g,a, b) around 
a = b = 0, and use symmetry arguments to dramatically 
reduce the number of independent terms at each order. 
When compared with published simulation results, this 
"spin expansion" seems to rapidly converge: the leading- 
order terms yield a surprisingly good first approximation, 
the ncxt-to-lcading-order terms give an even better ap- 
proximation, and so on. In the present paper, we present 
these points in detail, and explore some of their implica- 
tions and extensions. 



A. Some advantages of the spin expansion 

How does the spin expansion complement other ap- 
proaches to the BBH merger problem? 

1. First, it is simple. Previous approaches — notably 
the post-Newtonian approximation and full numerical 
relativity — are highly technical and sophisticated and 
have taken decades to develop. By contrast, we believe 
that the derivations in this paper will be accessible, even 
to physicists with no prior expertise in this area. 

2. Second, it is general. In this paper, we focus on 
applying the spin expansion to the final black hole's mass. 
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kick, and spin (m,k, s), but we expect that analogous 
arguments may be useful for studying any other final 
observable with well defined transformation properties 
under the simple symmetries {R, P, X} discussed below. 
This may also include the multipoles of the gravitational- 
wave signal emitted by the binary — this is currently a 
speculation, and remains a topic for future work. 

3. Third, it is conceptually distinct. Previous ap- 
proaches attempt to follow the dynamics of the merger — 
by solving the Einstein equations numerically, or through 
some analytical approximation. By contrast, our ap- 
proach is to consider the map directly from the initial 
state to the final state, and thereby "leap over" the com- 
plicated merger dynamics in between. To constrain this 
map, we rely purely on symmetry arguments, together 
with the assumption (supported by simulations) that the 
Taylor expansion of the map around a = b = con- 
verges rapidly. Therefore, our approach clarifies which 
aspects of the final state are due to the complicated non- 
linearities of Einstein's equations, and which aspects fol- 
low from more elementary considerations. 

4. Fourth, it is practical for cosmological and astro- 
physical applications. For example, a cosmological simu- 
lation of galaxy merger may also wish to track the corre- 
sponding supermassive black holes, since their feedback 
may be important for galactic structure and evolution. 
Of course, it would be hopeless to follow the BBH dy- 
namics in detail — the dynamical time near final merger 
is too short relative to the other timescales in the prob- 
lem. Instead, one is likely to resort to a simplifying algo- 
rithm: e.g. when the two holes get sufficiently close, they 
are replaced by a single hole with appropriate quantities 
{m, k, s}. The fact that the spin expansion maps the ini- 
tial state (well before merger) directly to the final state 
(well after merger) makes it well suited for these types of 
problems. 

5. Fifth, it is efficient. To fully specify the initial 
BBH spin configuration, we must specify 6 numbers — 3 
components each for the initial spins a and b. What is 
the most economical way to map out this 6-dimensional 
space with numerical simulations? If we crudely put 10 
grid points along each direction, we would need 10^ sim- 
ulations — an impossibly large value, since each simu- 
lation is very computationally expensive. On the other 
hand, we can use the spin expansion to map out the same 
6-dimensional space, at second or third-order accuracy, 
with only 0(10) simulations — a huge computational 
savings! This issue is treated in detail in Sec. |Vl In this 
sense, the spin expansion acts like a kind of "data com- 
pression," describing the 6-dimensional space of initial 
spins more succinctly, without oversimplifying it. 

6. Sixth, it is valid through the entire merger. The spin 
expansion is based on exact symmetries of general rela- 
tivity, which are equally valid during all stages of BBH 
merger: inspiral, plunge, and ringdown. Hence, it is par- 
ticularly well suited for questions about how the final 
state depends on the initial state in BBH merger. By 
contrast, the post-Newtonian approximation is an expan- 



sion of the Einstein equations in u/c, and treats BBHs as 
pairs of interacting point particles. It breaks down during 
the late stages of the merger — the last orbits, plunge, 
and ringdown, when the black holes begin to orbit with 
relativistic velocities and then cease to be two separate 
entities. Since these late stages emit gravitational waves 
copiously, and play a crucial role in determining the prop- 
erties of the final black hole, the post-Newtonian formal- 
ism is not well suited for predicting the final state of BBH 
merger. (On the other hand, it is excellently suited for 
describing the binary and its gravitational-wave emission 
when the black holes are well separated and orbiting non- 
relativistically.) 

7. Seventh, it is predictive. As we show in detail in 
this paper, the spin expansion makes a host of detailed 
(and successful) quantitative predictions for the results of 
simulations. Many of these predictions are new and dis- 
tinct from the predictions of other analytical approaches. 
They reveal interesting features of the simulations that 
might not have been noticed otherwise. 

8. These predictions are derived — they are not 
guesses. This is to be contrasted with expressions for 
the final kick velocity [1, i, i, ^ which were inspired by 
post-Newtonian equations, but are ultimately empirical 
fitting formulae which are not derived. Indeed, if these 
formulae continued to hold in general, it would be quite 
amazing. They are linear in the initial spins, and we will 
highlight several effects which appear to be quintessen- 
tially non-linear in the spins, and hence not captured by 
the post-Newtonian-inspired fitting formulae. 

In this subsection, we have made an aggressive case 
for the merits of the spin expansion. We must also stress 
the obvious point that the spin expansion merely comple- 
ments the other approaches to BBH merger — it does not 
replace them! It should be clear that numerical simula- 
tions, post-Newtonian techniques, and post-Newtonian- 
inspired fitting formulae offer a huge amount of dynam- 
ical information and insights which cannot be obtained 
from the spin expansion alone. Nonetheless, the spin ex- 
pansion provides unique understanding, as we hope will 
become clear in the course of this paper. 



B. Observational motivations 

The process of BBH merger is of great observational 
interest, since it is expected to govern the final evo- 
lution of both stellar-mass BBHs (produced in stellar 
collapse) and supermassive BBHs (produced whenever 
galaxies merge). 

Gravitational waves from merging stellar-mass BBHs 
are an important source for the ground-based detector 
LIGO [25| , while gravitational waves from merging super- 
massive BBHs are a primary source for the space-based 
detector LISA [2^. These gravitational waves will pro- 
vide unprecedented tests of strong-field general relativity, 
and open a new window onto exotic and previously invisi- 
ble astrophysical phenomena. For example, after the two 
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initial black holes form a common event horizon, they are 
predicted to "ring down" (like a bell) to a final quiescent 
Kerr black hole. This ring-down signal is an important 
observable for future gravitational wave detectors, and 
may be thought of as a superposition of so-called "quasi- 
normal modes." The observed spectrum of quasinormal 
modes depends on the quantities {m, k, s} characterizing 
the final black hole, so the ability to predict these quan- 
tities may play an important role in interpreting these 
observations. 

The quantities {m, k, s} characterizing the final black 
hole are also interesting astrophysically. These quantities 
may be probed observationally via the x-ray spectrum 
emitted from the inner edge of the accretion disk around 
a black hole. From the standpoint of supermassive BBH 
merger, they are believed to be linked to a variety of ob- 
servables, including: (i) the quasar luminosity function 
[27I . [28| : (ii) the location of a quasar with respect to its 
host galaxy [2^; (iii) the orientation and shape of jets in 
active galactic nuclei [13, HH; {iv) the correlation between 
black hole mass and velocity dispersion in the surround- 
ing stellar bulge [s^. [ssl; {v) the density profile in the 
centers of galaxies SJ, ■ Supermassive BBH merger is 
part of an interconnected web of astrophysical issues — 
the structure of this web is still poorly understood, and 
many fascinating questions remain. 

Finally, since the present paper is concerned with the 
merger of spinning black holes, we should note that many 
— perhaps even most — astrophysical black holes are 
indeed expected to have significant spin. On theoretical 
grounds, black holes grown by gas accretion are expected 
to have spins near the maximum Kerr limit [Sa . [sj- 
These predictions are supported by recent observations 
of active galactic nuclei (AGN) by the X-ray observatory 
XMM-Newton. Features in the Fe-Ka line of the Scyfcrt 
1.2 galaxy MCG — 06-30-15 are best modeled by a black 
hole with Kerr parameter a = 0.989to:oo2 at 90% con- 
fidence (where a = 1 would correspond to the maximal 
Kerr value) [3l|. Other observations suggesting signifi- 
cant black hole spin include [s^, [4^ . While torques from 
accreting gas tend to align the orbital and spin angular 
momenta with the large-scale gas flow in gas-rich "wet" 
mergers 14 1| . no such mechanisms exist in gas-poor "dry" 
mergers [43 making studies of generic initial BBH spin 
configurations essential for understanding these systems. 



C. Outline of this paper 



lations of these configurations and some less symmetric 
ones in Sec. IIVI Not only is the spin expansion consistent 
with all existing simulations, but it allows us to discover 
qualitatively new non-linear effects by specifying the spin 
dependence these effects must take. These non-linear 
spin effects reveal the inadequacy of existing "Kidder" 
fitting formulae for kick velocities 0, [l^l modeled after 
the purely linear terms in post-Newtonian expressions for 
the instantaneous loss of linear momentum. Though ex- 
isting simulations verify these exciting predictions of the 
spin expansion, they fail to fully constrain many of the 
terms at second order and beyond. In Sec.|Vl we propose 
a new series of simulations that will allow us to calibrate 
the coefficients of all terms to second order. We hope that 
the promise of our approach and its successes described 
in this paper will motivate numerical relativists to un- 
dertake these simulations in the near future. Finally, in 
Sec. IVIl we take stock of what has been accomplished in 
this paper and what remains to be done. 

Supplementary information has been organized into a 
series of appendices. We compile the relevant results 
of recently published simulations for convenience in Ap- 
pendix 1^ Lengthy equations relating different spin ex- 
pansions in Sec. IIVI arc relegated to Appendix iBl for clar- 
ity of presentation. Appendix [C] provides third-order 
terms in the spin expansion to extend the second-order 
calibration procedure described in Sec. El In AppendixiDl 
we show how systematic uncertainties in simulations af- 
fect our estimates of the coefficients in the spin expan- 
sion. In particular, for a general set of simulations of 
given covariance, we derive the minimum-variance unbi- 
ased estimators for the spin expansion coefficients. Fi- 
nally, in Appendix [El we briefly remark on the general- 
ization of the spin expansion to initially eccentric orbits. 



II. THE SPIN-EXPANSION FORMALISM 



For a brisk introduction to the spin expansion, see the 
first two pages of [l|. In this section, we introduce the 
same formalism at a more leisurely pace, providing more 
detailed explanations along the way.^ 



This paper is organized as follows. In Sec. [Hi we review 
the formalism of the spin expansion introduced in our 
Letter [l|. We use this formalism in Sec. IIIII to identify 
several particularly symmetric initial spin configurations 
for which a subset of the final quantities / G {m, s^, ki} 
identically vanish. These configurations may be useful to 
numerical relativists to help identify systematic errors in 
their codes that violate these symmetry constraints. We 
compare the predictions of our spin expansion to simu- 



If you find the presentation in this section suspiciously simplistic 
or Newtonian, we request your patience. Although we regard 
this simplicity as one of the key virtues of the spin expansion, 
in a forthcoming paper we make contact with the full 34-1 for- 
mulation of general relativity, and also with the post-Newtonian 
expansion. In the meantime, the proof is in the pudding: we 
hope that Section llVI will convince you that the spin expansion 
is powerfully explanatory and predictive when compared with 
existing simulations. 
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FIG. 1: A BBH system, including the orthonormal triad de- 
fined in the text. 



A. Preliminaries 

Imagine two black holes, A and B, in a circular orbit^, 
as shown in Fig.[T] Assume that A and B are initially far 
apart — far enough that they may be thought of as two 
Kerr black holes, characterized by masses {Ma, Mb) and 
spins (Sa, Sf,). Let us work in the center-of-momentum 
frame of the complete system (including the gravitational 
radiation, and the momentum that it carries). The or- 
bit gradually shrinks due to gravitational- wave emission, 
until A and B eventually merge. After the merger, the 
spacetime quickly settles down to a final Kerr black hole 
with mass M/, spin S/, and recoil velocity (or "kick ve- 
locity") k relative to the center-of-momentum frame. 

Next recall that classical general relativity has a trivial 
one-parameter rescaling symmetry. Starting from a given 
solution, we can obtain another solution by rescaling ev- 
ery physical quantity X according to its mass dimension 
dx'- X —> X^^ X , where A is an arbitrary positive num- 
ber.^ Of course, rescaling the solution in this way is 
closely related to rescaling the basic unit of mass. For 
our purposes, this freedom is more distracting than in- 
teresting: from now on, we will work exclusively with 
dimensionless quantities, which are unaffected by such a 
rescaling. In particular, it is useful to define the dimen- 
sionless initial mass ratio 

q = Mb/Ma, (1) 

the dimensionless initial spins 

a=SjM2 h=Sb/Ml (2) 



^ Gravitational radiation carries away energy more efficiently than 
angular momentum, and hence circularizes BBH orbits |43|| . 
Thus, most astrophysically relevant systems are expected to cir- 
cularize l ong before merger. With the exception of a few papers 
{e.g. [il, Ha, 0, 113, simulations of BBH mergers thus far 

have focused on circular orbits. Nevertheless, in Appendix [E] 
we will explain the extension of our formalism to non-circular 
(eccentric) orbits. 

^ Throughout this paper, we work in "geometrical units" with 
= c = 1, so that each physical quantity has units of mass to 
some power, called its "mass dimension." For example, angular 
momentum has mass dimension = 2, while distance and time 
both have mass dimension = 1. 



the dimensionless final mass 

m EE Mf/{Ma + Mb) (3) 
and the dimensionless final spin 

s^Sf/Mj. (4) 

B. Initial configuration of a blaclt hole binary 

To fully specify the initial configuration of this binary 
system, how much information do we need to provide? 

Since we are only interested in dimensionless quanti- 
ties, we can fully specify the initial state of the binary in 
terms of 8 numbers {V', 9, Oj, 6^} as follows. First choose 
an inspiral parameter ip, by which we mean a dimension- 
less quantity that varies monotonically along the orbit of 
the binary during the adiabatic inspiral phase. For ex- 
ample, "0 could be the (dimensionless) orbital separation 
r /{Ma+Mb), or the (dimensionless) magnitude of the or- 
bital angular momentum L / {Ma + Mb)'^ . At the "initial 
instant" (i.e. at some particular value of -0), define an or- 
thonormal triad {e*^^', e'^', e^'^^} as shown in Fig.[l] e'^' 
points along the orbital angular momentum, e'^' points 
from A to B, and e*^^) = e^'^) x e^^^. This same triad 
is conventionally introduced in post-Newtonian studies 
of spinning compact binaries (see e.g. [i^, [s^, HH 113] )■ 
Now we can specify the initial state of the binary, at the 
initial instant "0, by giving 7 more numbers — namely 
the dimensionless mass ratio q, and the initial spin com- 
ponents 

a, = a • e(^) = b • e^'^ (5) 

relative to the orthonormal triad. 

We can think of the 7 numbers {q,a^,bj} as parameter- 
izing the 7-dimensional space of physically-distinct black 
hole binaries (in circular orbit). Each "point" in this 7- 
dimcnsional space corresponds to an inspiral trajectory. 
On each trajectory, there is an "initial instant" labeled by 
■0, at which the 7 numbers {q, a^, b^} are to be specified. 

C. Symmetry considerations for binary systems 

Now let us consider how various final (post-merger) 
quantities can depend on the initial quantities presented 
in the previous subsection. 

Long after the merger, let / denote a final quantity 
of interest. For the purposes of illustration, we will fo- 
cus in this section on a particular set of final quantities 
/ £ {m, fcj, s^} — namely, the final Kerr black hole's di- 
mensionless mass m, and the components 

ki=k- e(') Si = s ■ e(*^ (6) 

of its final kick and spin — relative to the orthonor- 
mal triad defined at the initial time ijj as explained in 
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the previous subsection. Our presentation will hopefully 
be general enough to make it clear how to apply the 
same formalism to various other final quantities of in- 
terest, such as the (dimensionless) total radiated energy 
Ei-ad = ET-s.d/{Ma+Mb), or the (dimensionless) total radi- 
ated angular momentum jrad = 3^^^/{Ma + MbY . 

Any dimensionless final quantity / is a function of the 
initial quantities {q,a^,h,j} and the initial instant '0 s-t 
which these quantities were specified, 



(7) 



The goal of this section is to constrain the function 
9, Oj, 6j) as much as possible, using symmetry ar- 
guments alone. We will consider 3 simple transforma- 
tions of the binary system: rotation "i?," parity "P," 
and exchange "X." Once wc know how the initial and 
final quantities transform under i?, P, and X, we can 
constrain the map f{ip,q,a^,b^) by requiring it to relate 
initial quantities to final ones in a way consistent with 
their respective transformation laws. 

First consider the transformations R and P. i? is a 
global 3-dimensional rotation of the entire binary sys- 
tem (as if it were a single rigid body) , and P is a global 
parity transformation which reflects every point in the 
binary through the origin (the center of mass). How do 
the initial quantities {'0, q, a,, 6^} transform under R and 
PI The quantities {e'"'^', e'-^-'} are both vectors^ while 
the quantities {e(^',a, b} are all pseudovectors. The dot 
product of two vectors or two pseudovectors is a scalar, 
which is invariant under both R and P. The dot product 
of a vector and a pseudovector is a pseudo scalar, which is 
invariant under R and flips sign under P. Thus, the quan- 
tities {o]^, ^1, 62} are pseudoscalars, while the quanti- 
ties {03,63} are scalars. In summary, the initial spin 
components transform under P as 



where we have introduced the convenient notation 



hi = {-0^,-03, 
hi = {-61,-62, 



(8) 



(9) 



Note that the initial mass ratio q is also a scalar, and the 
parameter may be chosen to be a scalar: for example, 
the magnitude r/{Ma + Mb) of the initial separation, or 
the magnitude ^/(Afa+Mb)^ of the initial orbital angular 
momentum. 

For the rest of this paper, we restrict our atten- 
tion to final quantities / that are either scalars (like 
{?n, /cj^, fcjj S3}) or pseudoscalars (like {s]^, Sj, ^3}). In 
other words, we focus on final quantities / that are in- 
variant under i?, and transform under P as: 



./ - (±)p/, 



(10) 



where (±)p = +1 when / is a scalar, and (±)p = — 1 when 
/ is a pseudoscalar. The function f{tp, q, a^, 6^) automat- 
ically respects R (since the initial and final quantities are 



both invariant under R). In order to be consistent with 
P, it must satisfy the constraint 



/(0,g,a,,6,) = (±)p/('(/;,g,a,,6,). 



(11) 



Finally consider an "exchange transformation" X, 
which leaves the physical system absolutely unchanged, 
but simply swaps the labels of the two black holes, 
A B. How do the initial quantities {"0, q, a^, 6^} trans- 
form under X7 The parameter ip is invariant, and the 
mass ratio transforms as q ^ 1/q. The initial spin vec- 
tors arc swapped, a ^ b, while the triad elements trans- 
form as 



{eW,e(2),e(3)}^{-e 



(1) 



.(2) 



-e(3)} 



(12) 



Therefore the initial spin components transform as 



6,- 



(13) 



Now suppose that the final quantity / transforms under 
X as: 



(14) 



where (±)j^=-|-l when / is "even" under exchange (like 
{A:3,S3,to}), and (±)jc = — 1 when / is "odd" under ex- 
change (like {ki,k2, Si, 82}). In order to be consistent 
with X, the function f{tp,q,a^,b^) must satisfy the con- 
straint 



f{ijj,q,ai,bi) = {±)xfitp, 1/9, 6,, 5j). 



(15) 



Equivalently, but more conveniently, if / transforms un- 
der the combined transformation PX as: 



/ - i±)pxL 



(16) 



then the function f{ip,q,a^,b^) must satisfy the con- 
straint: 

f{ip,q,ai,bi) = {±)pxf{^IJ, l/q,bi,ai). (17) 
Note that (±)p^ is just given by the product 

{±)px = (±)p(±)x- (18) 

We emphasize that Eqs. (fTT|) and (fT7|) capture the key 
results of this subsection. To illustrate these formulae, 
let us apply them to the final quantities / g {to, k,^, s^}. 

First, in Table [H we collect the corresponding values 
of (±)p, (±)^, and {±)px for / G {m, fc,,sj. These 
values are easy to check. For example, consider the com- 
ponent Si = s ■ e^^\ Under a parity transformation P, 
the pseudovector quantity s (an angular momentum) is 
unchanged, while the vector quantity e*-^-* (the direction 
from A to B) flips sign, so their dot product Si is a pseu- 
doscalar: (±)p = —1. Under an exchange transformation 
X, which merely changes the labels A ^ B, the final an- 
gular momentum s is clearly unchanged, but the triad 
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/ |m, S3|si, S2|fci, fc2|fc3 

(±)p I + I I + h 

{±)x + z _± 

(±)PA- + + - - 



TABLE I: Transformation under P, X, and PX, for various 
final quantities /. 

element e'-^-' (the direction from A to B) flips sign, so 
their dot product is odd under exchange: (±)jf = — 1- 
Now, using Eq. (fTT|) . together with the (±)p row in 
Table m we find that parity P implies that the final quan- 
tities / € {m, k-, Sj} must obey the following constraints: 

to(V', q, flj, 6,;) ^ +m{ip, q, d.^, b.^) (19a) 

ki{ilJ,q,a„b,) = +kj^{il;,q,d,,b.^) 
fc2(V', 9, a,, 6J = +fc2(V', g, a,, &,) (19b) 
fc3(V', g, flj, ^j) = -fc3(V', 9, a,, &,) 
Si(V',g,a,;,6,;) = -.S;^(?/>,g,a,j,&,J 

S2{'fp,q,ai,K) = -S2(V',g,ai,&i) (19c) 
S3(V',9,a,;,6,;) = +s3(V',g,ai,&i) 

Using Eq. ([T7]), together with the (±)pjf row in Table 
m we find that exchange symmetry (or, more correctly, 
PX) implies that the final quantities / € {m, fc^, s^} must 
obey the following constraints: 

m{ip, q, a^,b^) = +m{ip, 1/q, 6,, aj (20a) 

ki{tp,q, ai,b^) = -fci(V', l/g,&,,a,) 

A:2(V', g, a,, 6.) = -^O', 1/g, a.) (20b) 

^3(V', 9, a^, ^j) = -fc3(V', 1/9, a^) 

Si(V^,(7,a.,6J = +Si(V', l/<7,fc„aj 

S2(V', (?, a,, 6,) = +S2(V', 1/9, a.) (20c) 

S3i'^,q,ai,b^) = +s^{ij,l/q,b^,a^) 

D. Series expansions for the final observables 

Symmetry considerations impose important con- 
straints on the maps f{^jj,q,a^,b^), but to make further 
progress it is useful to Taylor expand these maps about 
a = b = 0.^ In terms of the spin components {a^,b^}, 
this "spin expansion" can be written in the form 

/ = /™l"2™3l"l"2»3(^^q)a™la™2a^^6?;'l6^=^6^'^ (21) 



* In performing this Taylor expansion we assume that the map 
f{^p, q, o.^,b^) is analytic in the neighborhood of a = b = 0. While 
recent work by Pretorius and Khurana |44| suggests that certain 
finely tuned eccentric orbits can be exponentially sensitive to 
initial conditions, stable circular orbits of non-spinning BBHs 
should remain circular until the final plunge and merger. 



where separate summations over the 6 different indices 
{rU]^, mj, mg, ri]^, 7121 '^3} from to 00 are implied. Note 
that the expansion coefficients y™i™2™3|"in2n.'i g^.^ ^^^y^ 
independent of and 6j, but still depend on 1]: and q. 
Naive counting suggests that the number of terms in 
these expansions should grow rapidly with increasing or- 
der in the initial spins (1 zcroth-order term, 6 first-order 
terms, 21 second-order terms, and so on). However, the 
transformation requirements imposed by P and X signif- 
icantly reduce the number of terms that actually appear. 

Under a parity transformation P, the quantities 
{fl]^, flj, 6]^, 62} change sign, implying that individual 
terms in the expansion arc multiplied by (—1)''', where 
we have defined 

7 = nil + '^2 + '^1 + (22) 
Eq. (jlip therefore implies the Parity Constraint 

y™l™2™3l"l"2"3(q)=;(±)^(_l)7 j™im2™3l"l"2"3(q) 

(23) 

This may be restated as the Parity Rule: 

• Only terms with even (odd) 7 can appear in the spin 
expansion of a (pseudo) scalar f. 

Additionally, if we expand both sides of Eq. p7|) and 
equate terms with the same dependence on the initial 
spin components, we obtain the Exchange Constraint: 

J™l'"2™3l"l"2"3(^ — (-|-)^^,y"l"2"3l™l™2™3(^ X/g) 

(24) 

Let us again illustrate these results by applying them 
to the final quantities / £ {m, fcj, s,j}. We start by writing 
the expansion as 

m = m"'im,m,\n,n,n, (^^ ^)„mi „m2 „m3 ^„3 (25a) 

ki = fci™i™^'"=^i"i"^"=^(g, v)a"^arar^r^r^r 
^2 - fc2™'™''"''"'"'"'(9>'/')ar'arar^r^2'&r (^^t) 

fcg = fc3"i™="3l"i"^"3(g,V')a™^arar^r^2'^r 

Si = Si"i"='"3l"i"^"3(g,^)a™ia™=a™^61'i6^^&^'^ 

.S2 = S2'"''"''"''"'"'"'(9>V')ar'arar&r^2'&3' (25c) 
.S3 = S3™i"="=^l"i"^"=^(g, V)a™ia™^a™='65'i5^^&f 

The parity constraint, (fTT|) or implies that only 

terms with the correct parity (even 7) can appear in the 
expansions of the scalar quantities {m, k^, /cj, S3}. Terms 
with the wrong parity (odd 7) must vanish. Conversely, 
only terms with odd 7 have the correct parity to appear in 
the expansions of the pseudoscalar quantities {si, S2, fcg}; 
terms with wrong parity (even 7) must vanish. 

The exchange constraint, ([TT]) or implies that 

the remaining coefficients must satisfy the following con- 
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straints: 









™l™2™3(l/g) 


(26a) 




"l"2"3 {^q) = 


^ "l"2"3 


™1™2™3 




"-2 


"l"2"3 ^g^„ 


, ^1*^2^3 

"-2 


(1/9) 


(26b) 




"l"2"3(^^-)^ 


, "l"2"3 


(1/9) 




^1 


"l"2"3(-gj_ 




™1™2™3 




^2 


n-, n^n.:. / \ 

(9) = 


"l"2"3 


'(1/9) 


(26c) 


^3 


"l"2"3J^^j_ 


+ S3 


'™1™2™3 





where, for brevity, we have not displayed the ip depen- 
dence on each side. 

Note that the exchange constraint imposes dual- 
ity relations between coefficients at q and 1/q. Without 
loss of generality, we can focus on the region of initial 
parameter space with < g < 1, since spin expansions in 
the region with 1 < g < c» may simply be obtained via 
the duahty relations or p5)) . 

In the equal-mass case, q = 1 = 1/q, so the 
duality relations directly relate previously indepen- 
dent coefficients. In particular, these relations re- 
quire the kick- velocity coefficients ;j.™i™2™3l"i"2"3 ^f^i^]^ 
{mi, 777.2, 7773} = {771, 772, '^3} to vauish for q = 1 since the 
components fcj have (±)px = — 1- This result is consis- 
tent with the famous kick formula of Fitchett [s^, who 
was one of the first to calculate the gravitational-wave 
kick resulting from the merger of non-spinning BBHs in 
the Newtonian approximation. In his honor, we would 
like to name the non-spinning kick coefficients f^^^^^^^^ 
"the cocFitchetts." 



III. EXACT RESULTS AND SPECIAL 
CONFIGURATIONS 

In this section, we will highlight several exact results 
which follow from the symmetry considerations devel- 
oped in the previous section. Since they arc exact, these 
predictions may be useful for testing numerical codes 
which compute BBH mergers numerically. 

To derive and summarize the results, it is helpful to 
think about the operations P and X as elements of a 
discrete group G of operators acting on binary systems S. 
This approach conveniently generalizes to include other 
discrete symmetries like charge conjugation C.^'^ The 





P 


X 


PX 


c 


PC 


xc 


PXC 


q 


q 


1/9 




q 


q 






ai 






+K 


+ffli 


-0.1 


-bi 


+bi 


a-2 


-02 


-h2 


+b2 


+02 


—02 


-b2 


+&2 


"3 




+&3 


+63 


+03 


+0.3 


+63 


+63 


&i 




-ttl 


+ai 


+bi 


-bi 




+ai 


&2 


-b2 


-^2 


+an 


+62 


-b2 


—02 


+0.2 




+ &3 




+ 0.3 


+b3 


+b3 


+03 


+^3 


Qa 


+ Qa 


+Qb 


+ Qb 


-Qa 


-Qa 


-Qb 


-Qb 


Qt 


+ Qb 


+ Qa 


+ Qa 


-Qb 


-Qb 


-Qa 


-Qa 


m 


+777 


+ 777 


+ 777 


+777 


+ 777 


+77t 


+777 


ki 


+ fcl 


-fcl 


-fcl 


+ fcl 


+ki 


-fcl 


-ki 


^2 


+k2 


— ^2 


-k2 


+ k2 


+k2 


-k2 


— ^2 


^3 




+^3 


-k3 


+k3 


-k3 


+k3 


-k3 


Si 




-Si 


+ Sl 


+Sl 


-Si 


-Si 


+Sl 


S2 


-S2 


-S2 


+ S2 


+S2 


-S2 


-S2 


+ S2 


S3 


+ S3 


+ S3 


+ S3 


+S3 


+S3 


+S3 


+S3 


Qf 


+Qf 


+Qf 


+Qf 


-Q.f 




-Q, 





TABLE 11: Transformations of the initial and final observables 
listed in the first column under the group of operations G 
formed from parity P, exchange X, and charge conjugation 
C and listed in the first row. 



three operators P, X, and C all square to unity and 
commute: 



P2 =X2 =C2 = 1, 

[P,X] = [P,C] = [X,C]=0, 



(27) 



so they generate the Abelian group G = Z2 x Z2 x Z2, 
with eight elements:^ 



G={l,P,X, PX, C, PC, XC, PXC}. 



(28) 



The non-trivial elements of this group transform the ini- 
tial and final quantities in a BBH merger as summarized 
in Table HH 

Each column in Table |TT] (beyond the first) is identified 
with an operator g G G, and establishes a relationship 
between two physically distinct BBH systems S and S' 
related by S' = gS, S = gS' . For example, the second 
column (P) relates any system S with initial spin com- 
ponents {flj, 6j} to a second system S' with initial spins 



{a'^,a'^,a'^} = {- 



{b\,h'^,b'^}^{-b„-h^,+h^] 



(29) 



Charge conjugation C is an exact symmetry of classical general 
relativity and elcctromagnctism. Under C, the charge Q of a 
black hole changes sign, Q ^ ~Q. 

Though black hole charge is expected to be negligible in astro- 
physical contexts (and certainly much smaller than the Kerr- 
Newman bound) , it could be important in microphysical contexts 
— for example if TeV-scale quantum gravity leads to black-hole 
production in high-energy accelerators, such as the Larg e Hadron 
Collider (LHC), soon to begin operation at CERN [5411 . 



^ Time reversal T is another exact symmetry of classical general 
relativity and elcctromagnctism. One might hope that we could 
use T to impose further constraints on binary black hole merger, 
but unfortunately we cannot. Black hole merger is inherently 
a dissipative process; we are interested in black holes that emit 
gravitational waves and inspiral — not those which absorb grav- 
itational waves and outspiral! 
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input 


output 


p 


a-, = = b-i = 6t = 


Si = So = it-j = 

■^1 "^2 ""3 ^ 


X 


(a-, ttol = — (b^ bn) 

a3 = bs, q=l, Qa = Qb 




PX 


(n.^ n,n n.o^ = (h^ hn ho^ 
; ^2 ; "'3/ ) i 

g = 1, Qa = Qb 


= /j2 — ^3 — 


c 




o,=o 


PC 


ai=a2 = &i = &2 = 0, 


Sl=S2 = fc3 = Q/=0 


xc 


(ai.aa) = -(&i,&2), 
a3 = &3, g=i, Qa = -Qb 


fci=fc2 = Si = S2 = Q/ = 


PXC 


{ai,a2,a3) = (6i,b2,fe3), 

g=i, Qa = -Qt, 


fcl = ^2 = ^3 = Q J = 



TABLE III: Special initial configurations ("input"), and the 
corresponding predictions for tlie final state ("output"). 



The final quantities for the system S' will be given in 
terms of those for S by 



{m',Q'^.}^{m,Qf.}. 



(30) 



Some of the 7 predicted relationships represented by the 7 
columns of Table Hi] may be useful to numerical relativists 
for identifying errors in their numerical codes that fail to 
respect the given symmetries. For example, P demands 
that simulations of the systems S and S' must yield final 
quantities that are related by Eq. ([30|l . Any deviations 
from this relationship would reveal systematic errors in 
the simulations that need to be corrected. 

Each non-trivial element .g G G not only establishes a 
dual system S' = gS for all BBH systems 5, but also de- 
fines a "special" configuration Sg that is its own dual. i.e. 
gSg = Sg. Since the initial state is invariant under the 
final state must also be invariant, so we can conclude that 
any final quantity / that is not invariant under g must 
vanish. For example, consider the second column (P) in 
Table [III We see that, if the initial configuration satisfies 
^1 ~ ^2 ^ ^1 ^ ^2 ~ 0, then a parity transformation 
P leaves this initial configuration invariant. Therefore, 
we can conclude that the corresponding final state must 



satisfy = S2 



0, since these 3 quantities are not 



invariant under P. Similarly, from each of the 7 columns 
in Table |TT1 we can read off a "special" initial configura- 
tion, and the corresponding predictions for the final state 
of the system. These 7 special configurations, and their 
consequences, are summarized in Table Hill 



IV. TESTING OUR EXPANSIONS WITH 
EXISTING SIMULATIONS 

Although the previous section's exact results are in- 
teresting, the real power of the spin expansion is in the 



much larger set of approximate predictions it makes for 
generic spin configurations. 

Currently, published simulations of BBH mergers can 
be subdivided into 5 different classes of initial spin con- 
figurations. For each of these classes, we compare in de- 
tail the predictions of the spin expansion with existing 
numerical results. As we shall see, the spin expansion 
makes new and successful predictions in each case, and 
provides a simple and systematic way of deriving features 
of existing simulations that were previously opaque. 

Since rigorous systematic errors are not yet available 
for many of the simulated data sets analyzed in this sec- 
tion, our analysis of these simulations will by necessity be 
correspondingly non-rigorous. In particular, we will use 
the following rather heuristic procedure. For each data 
set, we compute the pcr degree of freedom (x^/d.o.f.). 
We do not attribute meaning to the overall value of the 
X^/d.o.f. (since, in many cases, we have had to guess er- 
ror bars for the simulated data). However, if including a 
new term predicted by the spin expansion leads to a large 
fractional reduction in the x^/d.o.f., we interpret this re- 
duction as admittedly non-rigorous evidence for this new 
term. Note that an overall rescaling of the error bars will 
not lead to a fractional reduction in the x^/d.o.f., so the 
conclusions in this section should be largely insensitive 
to our estimates for the error bars. 

In any case, statistical rigor is not the point of this 
section. Our goal is to illustrate our formalism through 
a few worked examples, to demonstrate its power to ex- 
plain currently available simulations, and to make several 
predictions for future simulations. It will usually be clear 
(by eye) that our leading-order and next-to-leading-order 
formulae provide a good explanation of the basic quali- 
tative features seen in the simulated data thus far. 

We have summarized our use of simulations in Ap- 
pendix [XI where we also explain our estimates for the 
corresponding error bars - which are often just crude 
guesses. The crudeness of these guesses sometimes leads 
to impossibly small x^/d.o.f. for our fits. This is pri- 
marily because genuine errors in the simulations are sys- 
tematic, whereas published works (and our own analysis) 
treat these errors as statistical. Systematic errors that 
preserve the symmetries of the configuration will be well 
fit by our expansions, regardless of their size, albeit with 
erroneous numerical values for the best-fit coefficients. 

In this section, we will often encounter 3-component 
quantities {xi,X2,x^) where we wish to treat the "1" and 
"2" components together, without the "3" component. 
Thus, it is convenient to introduce the notation 



{Xi,X2) 



(31) 



as these components are perpendicular to the orbital an- 
gular momentum vector. So, for example, aj^ = {a-^^,a2), 

= (/fc™^l°°°,fcr'™°)- From 
acts like a "2- vector" in the 



001 000 



kj^ = (ki,k2), and k 
a notational standpoint, 
sense that 



y_L =2^1 yi+ 2:2^2 ■ 



(32) 
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A. Case #1: Non-precessing spins, 

{q = 1 and a oc b oc e''^') 





First consider the case in which the black holes have 
equal mass {q = 1), and both spins are aligned (or anti- 
aligned) with the orbital angular momentum: 



(33) 



This corresponds to the special configuration identified 
with the operator P in the previous section, for which 
fcg and s_|_ exactly vanish. What about Tion-vanishing 
observables like kj_, S3, and m? 

The zeroth-order terms in the S3 and m expansions 
are non-vanishing, and physically correspond to the spin 
and mass of the final black hole produced in the merger 
of two non-spinning holes A and B. However, exchange 
X requires the zeroth order term in the kj^ expansion 

(the "coeFitchett" 1^™"'°"°) to vanish, as discussed in 
section Hi Dl To capture the leading-order (LO), next-to- 
leading-order (NLO), and next-to- next-to-leading-order 
(NNLO) behavior, we therefore expand Sg and m to sec- 
ond order in the initial spins and kj^ to third order. 

Using the parity and exchange constraints, (j23p and 
p4)) . to equate or eliminate coefiicients, our expansions 
([25]) for this configuration become: 



, , 001|000. u \ , 1 002|000/ 2 u'2\ 



+ kr'™°(ai-6i)-fkr'"^(ai63-6la3), 



000|000 



001 000/ , , X 

-■S3 [as + bs) 



002|000/ 2 , 72\ I 001|001 , 

•S3 (a3 + 03) + -S3 0303, 

^0001000 _j.^001|000(^^^^^) 
^002|000(„2^^2)^^001|001^^^^_ 



(34a) 
(34b) 
(34c) 



Unfortunately for us, currently published simulations 
only report the magnitude |k^| = kf + A:| , not the 
individual components ki and Taking the magnitude 
of the expansion (|34aP for k_|_ and Taylor expanding to 
third order in the initial spins yields 

1, I _ 1, OOIIOOOI I _, I 

X l+A{a.^ + b^) + B{al + bl) + Ca3b3 . (35) 

For convenience, we have introduced new coefficients A, 
B, and C, which may be expressed in terms of the origi- 
nal expansion coefficients ]j™i™2m3|nin2ii3^ These expres- 
sions are given in Appendix IB II 

It is interesting to note that, although |kj_| is even un- 
der both P and X (just like S3 and m), its expansion (|35p 
is different from the S3 and m expansions (|34bl[Mc|) . This 
is a reflection of the fact that, although |kj^| has the same 
transformation properties as Sg and m, it is a composite 



quantity constructed from more "fundamental" quanti- 
ties (k^ and k2) with different transformation properties. 
Looking more closely, we notice that the expression inside 
square brackets in Eq. ([55)) b close formal resem- 

blance to the expansions ((34bl[34cl) : there is a constant 
term at LO, a term proportional to (og + ^g) at NLO, and 
two terms proportional to (03 + bf) and a^b^ at NNLO. 
The compositeness of |kj^| manifests itself through the 
overall factor of |ag — 63 1 in front of the square brackets 
13. 



Eqs. ([35]) . (|34bp . and (|34cp make detailed quantititive 
predictions for the final kicks, spins, and masses — let's 
see how they stack up against actual simulations! 

Many groups have published results from the simu- 
lations of binary mergers with initial spins aligned or 
anti-aligned with the angular momentum direction e*-^-* 
d, [H, M, [H, [13, Hi- While all groups begin their sim- 
ulations with the binary on a quasi-circular orbit, they 
make different choices for the initial dimensionless orbital 
separation r/(M^ + Mf^). As explained in Section |TT1 our 
expansion coefficients (|25p are defined at a fixed value 
of the inspiral parameter "0, which in this case is the 
dimensionless orbital separation r/{M^ + M^). Connect- 
ing simulations performed at different values of "0 will in 
general require careful use of post-Newtonian approxima- 
tions as discussed in a forthcoming paper. However, for 
the special case considered here (where a cx b oc e^'^^), 
the initial spins will not precess and their projection onto 
the orthonormal triad {6*^^^ e*^^\ e'^'^^} will not vary with 
orbital phase. As such, it is possible in principle to 
jointly fit all the simulations even though they do not 
all correspond to the same initial separation. In prac- 
tice, there are systematic differences between numerical 
codes * which make a combined fit to all existing simula- 
tions unreliable. In Appendix I A 11 we describe our choice 
of simulations to test Eqs. ([55l I34bl I34cp . 

First consider the final kick |k, |. At leading order. 



Eq. (|35|) predicts that |k_|_| should be proportional to 
I13— 63I; this approximately linear behavior has been no- 
ticed in simulations in several previous papers [l^ . 
At next-to-leading order, Eq. ([55]) predicts that |k^| 
should receive a small additive correction proportional 
to 1 03 — 631(03 -I-63). This prediction is well supported by 
the simulations in the following sense. When we fit the 
28 simulations in [l^ with non-zero values for |k^ | to the 
leading-order term in Eq. ([35)) . there is only one fitting 
parameter (namely the magnitude |k']'^^'''''''|), and the fit 
is rather poor xVd-o.f. = 38.2/(28 - 1) « 1.4. Then, 
when we include the next-to-leading-order term, there is 
one more fitting parameter (namely A), and the fit im- 
proves dramatically: x^/d.o.f. = 5.14/(28-2) w 0.2. We 



For example, PoUney et al. Il4l discusses the necessity of prop- 
erly choosing an integration constant corresponding to the linear 
momentum acquired before the simulations begins. RezzoUa et 



ble for the discrepancy between their results and those of 



pon s 

m 
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Data set 


Best-fit coefficients 




|kxl 


H^001|000| 


PS 221 km/s 




(28 simulations) 


A « -0.205 






B f» -0.091 


m 




C ^ -0.201 


km 


S3 


0001000 
S3 R 


i 0.6893 




(48 simulations) 


0011000 
S3 ^ 


i 0.1524 






0011001 
S3 R 


i -0.0195 






0021000 
S3 R 


i -0.0121 




m 


^000(000 


=^ 0.9530 




(24 simulations) 


^0011000 


=a -0.0167 






^ooiiooi 


^ -0.0083 






^0021000 


=i -0.0052 





TABLE IV: Case #1 best-fit parameters, from fitting 
Eqs. pSI I34bl I34cp to the final kicks, spins, and masses of 
the simulations described in Appendix lA II These fits are 
displayed in Fig. [51 



interpret this significant drop in the x^/d.o.f. as strong 
evidence for the second-order term in Eq. ((35)) .^ This 



NLO fit, with only two fitting parameters (Ik^l and A), 
is displayed in the top panel of Fig. O 

Is there also evidence for the next-to-next-to-leading- 
order (NNLO) terms in Eq. ((35)) ? When we include 
these final two terms, there are two more fitting parame- 
ters {B and C), and the fit again improves significantly: 
xVd.o.f. = 2.2/(28 - 4) w 0.09. Our best-fit values for 
the 4 fitting parameters in Eq. (|35p are shown in Ta- 
ble ITVl 



Next consider the final spin Sg and the final mass m. 
It is convenient to treat these two quantities in parallel, 
since their expansion (j34bp and (|34cp are formally iden- 
tical to each other. Eqs. (j34bp and (j34c[) predict that at 
zcroth order the final spin S3 and final mass m should 
be equal, respectively, to the final spin g!^"'^!™'' and fi- 
nal mass TTjOOOIooo from the merger of two non-spinning 
black holes. At first order, there should be a linear cor- 
rection proportional to (03-1-63), and at second order there 
should be two corrections: one proportional to 0363, and 
the other proportional to {a'^ + b'^)}'^ 



PoWney et al. O 

also noted a significant deviation from linearity 
in the kick magnitudes for initially (anti-)aligned spins, although 
their guess for the correction term differs from that derived via 
our formalism. RezzoUa et al. [TBI arrived at the same second- 
order fitting formula as we did using similar considerations; our 
parameters jj-'JOilOOOj ^ correspond to |ci| and C2/C1 in their 
notation. 

RezzoUa et al. [Tsl argued that BBHs with equal and opposite 
spins (a3 = —63) should behave as if they were non-spinning, 
so that only a single term P2{o,^+b^)^ should appear at second 
order. According to this conjecture our coefHcicnts should be 
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FIG. 2: Case #1 best-fit curves. The top, center, and bottom 
panels show, respectively, the kick velocity |kx|, dimension- 
less spin S3, and dimensionless mass m of the final black hole. 
The data points appearing in each panel are described in Ap- 
pendix|Xl] The curves in the top, center, and bottom panels, 
correspond, respectively, to Eqs. (|35|) . ()34b|l . and p4c|l . with 
the coefficients given in Table IIVI Each curve has a fixed 
value of bs given by the legend at the bottom of the figure, 
with as varying along the abscissa. For presentation purposes, 
we have switched as and 63 for points on the magenta (long- 
dashed) curve; this exchange does not affect the values of |kx | 
and S3. 



These predictions are again supported by the simula- 
tions. First consider S3. By itself, the leading-order term 



3), the fit improves 
- 2) = 0.366]. Fi- 

002|000^ 2 , u2\ 1 

3 {ai+bi) and 



related as s„ 



1/2 



- P2- 



mo\om gj^^g ^ p^^^. [x^/d.o.f. = 3,661/(48 - 1) = 
77.9] to the 48 simulations of [11] and [13. When we 
include the NLO term, ^^^''^'^''(aa + b 
dramatically [xVd-O.f. = 16.8/(48 
nally, when we add the NNLO terms, s 
^ooi|ooi^^^^^ the fit is even better [x^/d.o.f. ~ 0.0386]. 
The fit to m using Eq. (|34cp is closely analogous: the 
zeroth-order fit is again lousy (x^/d.o.f. = 436/(24—1) = 
19.0); the first-order fit is much improved (x^/d.o.f. = 
46.8/(24 - 2) = 2.13); and the second-order fit is bet- 
ter stiU (xVd.o.f. = 5.30/(24 - 4) = 0.265). Thus, in 
both the S3 and m data sets, we have clear evidence for 
both next-to-leading-order (NLO) and next-to-next-to- 
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leading-order corrections.^^ Our best-fit parameters are 
shown in Table IIVI and the corresponding fits to Sg and 
TO are plotted in Fig. [21 in the middle and bottom panels, 
respectively. 

Though we have used numerical simulations to cali- 
brate the values of the coefficients in our spin expansions, 
physical intuition provides some insight into these val- 
ues. Test particles with orbital angular mometum aligned 
with the spin of a Kerr black hole have innermost stable 
circular orbits (ISCOs) with smaller radii, and therefore 
emit more gravitational radiation before merger. Similar 
behavior has been observed for comparable-mass black 
holes in numerical simulations [l^. This suggests, as 
we have indeed found by comparing to simulations, that 
the coefficient 77),ooi|ooo should be negative: more orbits 
implies more energy carried away in gravitational waves 
and a smaller final mass. The impact on the final spin is 
more ambiguous; aligned BBHs should convey more spin 
angular momentum to the final black hole, but less or- 
bital angular mometum because of their smaller ISCOs. 
Since Sg''^''^'"' is positive we conclude that the former ef- 
fect dominates over the latter though its comparatively 
small value evidences significant cancellation. 



expansion takes the form 

oo i 

/ = E E sin(j0) + 7('-^)a' cosO'^)] , (37a) 

where / represents one of the non- vanishing final observ- 
ables {e.g. fcg, Sg, or to). The coefficients ^/^''J) and 
■^Z^*'-''' are finite linear combinations of the original spin- 
expansion coefficients j™i™2m3|niTi2n3 ^ These linear com- 
binations can be derived explicitly by equating the two 
expansions (j25p and (|37a[) term by term; the first few 
combinations are provided in Appendix IB 21 

For situations in which varies while a remains fixed, 
we should go one step beyond Eq. (|37ap by collecting all 
terms with the same (j) dependence. Then Eq. (j37ap is 
rewritten in the form 

oo 

/ = E ['/^'^ sin(j0) + 7(^'^ cos(j(/.)] , (37b) 



B. Case #2: "Superkick" configuration 

(g = 1, a = -b, as = &3 = 0) 

We next consider the case in which the black holes have 
equal mass (g = 1) and equal and opposite spins lying in 
the orbital plane: 








(36) 



As shown in Section HITl (Table IIIip . this is an example 
of the special configuration associated with the exchange 
operator X, for which and Sj^ exactly vanish. What 
about the non- vanishing quantities /cg, Sg, and to? 

The "superkick" initial spin configuration ([36)) is pa- 
rameterized by only two numbers (a and (p), instead of 
six (flj and b^). As a result, when we substitute Eq. (j36p 
into the original spin expansions of Eq. (j25[) , many of the 
terms become degenerate with one another. It is conve- 
nient to eliminate this degeneracy by collecting all terms 
with the same dependence on a and (/). Then the spin 



where we have defined the coefficients 



3)r 



fU) 



E^/ 

i=3 



(38) 



We stress that Eq. (|37ap is the same series as Eq. (j37bp . 
The first form p7ap is appropriate for situations where 
a and (f> both vary, whereas the second form ([37bp is 
appropriate for situations where cj) varies but a does not. 

Symmetry considerations further restrict the terms 
that can appear in the expansions (j37p . Applying 
parity P to the initial configuration (j36p is equiva- 
lent to the transformation ch (h + n which sends 



{sin(j0),cos(j(/))} 
follows that 



{(-l)Jsin(j0),(-l)J"cos(j(/))}. It 



• In the superkick configuration iS6]) . only terms with 
even (odd) i and j can appear in the expansions 
|,?7|) for a scalar (pseudoscalar) quantity f . 



The final spins and masses are equally well fit by a single second- 
order term proportional to (ag-l-fcg)^ as suggested by RezzoUa et 
al. |15| . Though post-Newtonian approximations may suggest 
that the spin-spin term proportional to a^b.^ comes in at higher 
order, we restrict ourselves to what can be assorted purely on 
the basis of our formalism. 



Since the superkick initial spin configuration (j36p is in- 
variant under X, an exchange transformation does not 
yield any further constraints. 

These simple considerations lead to detailed quantita- 
tive predictions. To illustrate this point, let us start by 
displaying some of the leading terms in the expansions 



12 



(j37ap for the non- vanishing obscrvablcs fcg, Sg, and 



i-^^'^^fli + '^fcf '^'03^0(05)] sin (0) 
ci.(i.i).i_Lct,(3a)„3, 



' fcf '^^ _^ fcf flS + O (a^ )] sin (30) 



ci,(3,3)„3,c, (5,3)^5 



+ "4 ■^^a5+0(a^)]cos(3<?!)) 



(39a) 



rc„ (0,0)^0 ,cf2,0)^2 



+ 0{a^)] cos(O0) 



s„(4,2)^4 



a^ + Oia")] sin(2(/)) 

+ rSg^^'^^a^ ^ c^_(4,2)^4 ^ 0(^6)] cos(20) 



(39b) 



m= [^m("'°)a° + "m(2^°)a2 + 0(a^)] cos(O0) 

+ ... (39c) 

In these equations, we have exphcitly displayed the 
cos(0(/)) = 1 factors to emphasize the similarity between 
the j = and j ^ terms. 

First consider the predictions of Eq. (|39ap for the kick 
velocity fcg. The leading-order {0{a^)) terms predict that 
this kick should vary as sin(0 + phase) at a fixed value 
of a, and that the amplitude of this sinusoid should scale 
linearly with a. The ncxt-to-lcading-order {0{a^)) terms 
predict that this leading sin(0 + phase) behavior should 
receive a small additive correction of the form sin(3(/) + 
phase) with amplitude proportional to a^ . Next consider 
the predictions of Eq. (j39bp for the final spin S3. The 
leading-order (©(0°)) term predicts that the final spin is 
a constant, independent of both a and </>: S3 oc a° cos(O0). 
The ncxt-to-lcading-order (©(a^)) terms predict that the 
leading cos (00) behavior should receive a small additive 
correction of the form sin(2(/) + phase) with amplitude 
proportional to a^. Finally, since Eq. p9cp is formally 
identical to Eq. (j39bp . m should behave in the same way 
as Sg . 

We test these predictions using the published simula- 
tions of @, 01 1 with relevant data and errors described in 
Appendix lA 21 Although both papers estimate the final 
kicks fcg , the simulations in Q start at a different orbital 
separation r/{Ma + Mb) from those in 0. In Sec. lIVAl 
we were able to jointly fit simulations with different ini- 
tial separations because the relevant spin-expansion coef- 
ficients were insensitive to the initial separation. Unfor- 
tunately, in the present configuration p6p . the relevant 
expansion cocfhcients are sensitive to the initial sepa- 
ration. As we discuss in a future paper, it should be 
possible to connect the expansion coefficients at different 
initial separations using post-Newtonian techniques. For 
the time being, though, we must perform separate fits for 
the simulation sets {Ai} of and {Bi} of Q. As the 
initial spin magnitude a is fixed within each data set, we 
should use the spin expansion in the form (j37bp . Thus, 
for the non- vanishing obscrvablcs fcg, Sg. and m, we have 



Data Set 


A 


A 


B 


B 


^kl 


-350 


-323 


3753 


2714 


'kl 


1876 


1837 


-343 


-246 


^kl 


— 


-29.4 


— 


-20.8 






2.73 




83.5 


%jr.i/My 


0.2471 


0.2471 


— 


— 




— 


-0.0012 


— 


— 


"(Jrad/M^)' 




-0.0027 








— 


— 


0.6895 


0.6895 




— 


— 


— 


-0.0038 


a „2 








-0.0021 


'^(%-Brad)'' 


3.587 


3.600 






'^(%-Erad)^ 




-0.0301 






''(%-Erad)^ 




-0.0695 







TABLE V: Fits for Case #2. The first column lists the co- 
efficients being determined. Kick velocities are in units of 
km/s while the remaining coefficients are dimensionless. The 
second column provides the best-fit values for these coeffi- 
cients when the lowest-order terms are fit to data set A, the 
simulations of Expansions for the radiated angular mo- 
mentum Jrad /M'^ and energy %i5rad have zeroth-order terms 
while first order is lowest for k^. The third column lists best- 
fit values for next-to-lowest order fits; this is second order for 
Jrad/M^ and %-E'rad and third order for ^3. The fourth and 
fifth columns list the corresponsing values of coefficients for 
data set B, the simulations of [J. This data set provides the 
spin S3 of the final black hole instead of Jrad/M^ and %-Erad. 
The expansion for S3 has zeroth and second-order terms. 



the expressions: 

k^ = ''kl^hm{(j)j^ "-3 

-h'/fcf ^ sin(30) + '=A:f ^ cos(30) + 0{a^), 



"fc^^^ cos(0) 



(40a) 



S3 = '=S3(°^ + 's^^"' sin(2(/))+ "sg^"^ cos(20) + O(a^), (40b) 
TO = =TO(°'-h"m(2)sin(20) + '=m(2)cos(20)-he'(a'') (40c) 

Instead of reporting Sg and m, Campanelli et al. re- 
port the percentage of the initial energy carried away by 
gravitational radiation, %-E'rad, and the dimensionless ra- 
diated angular momentum J-^g^jM"^. Since both of these 
quantities are scalars, with P = -1-1 and X = -1-1, their 
expansions are exactly analogous to the expansions for m 
and Sg in Eqs. (|40b[) and (|40cp . The best-fit coefficients 
from our leading-order and next-to-leading-order fits to 
these quantities are listed in Table |Vl and the leading- 
order fits themselves are displayed in Fig. [3] 

First focus on the kick velocity fcg. The top panel in 
Fig. [3] clearly reveals the sin(0 -I- phase) behavior pre- 
dicted (at leading order) by the spin expansion, and pre- 
viously noted in Eq. (1) of d] and Eq. (6) of 0. Al- 
though the spin expansion cannot predict the phases of 
the sine waves in this panel, it predicts that (at lead- 
ing order) their amplitudes should be proportional to 



(2 



(2), 
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FIG. 3; Case #2 best-fit curves. Panels from top to bottom 
show tlie kick velocity ks, the radiated angular momentum 
^rad, the final spin S3, and the percentage of radiated energy 
-Erad, all plotted against the angle ((> between the initial spin 
a and e'^-*. The blue curves show fits to the square data 
points taken from 0], while the red curves show fits to the 
triangle data points of 0. The curves for fcs only show the 
first-order terms in Eq. (|4Uap . while those for Jiad, S3, and 
i5rad include all terms up to second-order. The x^/d.o.f. and 
best-fit coefficients are listed in Table IVl 

the spin magnitude a. This prediction is also supported 
by the simulations: the ratio of the best-fit values of 

[(^4^^)^ + {""k^^^??^^ for the simulations {Ai} of Q and 
{Bi} of is 0.684, not far off from the ratio of their 
spins 0.515/0.723 = 0.712. 

The next-to-leading order predictions of the spin ex- 
pansion for fcg also appear to be confirmed. To see 
this, we subtract the leading-order sin((/) -I- phase) pre- 
diction from the simulated data, and plot the resid- 
uals in Fig. m As predicted, these residuals oscillate as 
sin(3(/)+phase). Furthermore, the amplitude of this resid- 
ual oscillation scales as a^, as predicted. The ratio of the 
best-fit values of [{"k^^'')^ + ("fcf ^)2]i/2 for the simula- 
tions {Ai} and {Bi} is 0.343, which is close to the spin 
ratio cubed, (0.515/0.723)^ = 0.361. 

Thus, we have seen that the spin expansion succeeds in 
reproducing previously observed aspects of the k^ data, 
and also in predicting previously unrecognized features. 

The remaining three panels of Fig. [3] show observables 
whose 0-dependence agrees closely with what we expect 
for scalars like S3 and m in Eqs. (|40bp and (|40cp . At 
zcroth-order they are independent of 0, while at next-to- 
leading order the predicted sin(20 -I- phase) contribution 



100 



B 




-100 I ' ' ' 1 ' ' ' ' ' ' ' L. 

2 4 6 

^(radians) 



FIG. 4: The residuals Aks after the first-order terms of 
Eq. (|40a|l are subtracted from the simulated final kicks. As 
in the top panel of Fig. O the blue curves show fits to the 
square data points taken from 0], while the red curves show 
fits to the triangle data points of 0. These curves consist of 
the third-order terms of Eq. (|40a|l with the coefficients listed 
in Table lYl 



appears. Unfortunately for us, Q and chose to pub- 
lish different observables so we cannot test whether the 
amplitude of these double-frequency terms really scales 
like as predicted. Hopefully, future simulations per- 
formed at different values of a will address this question. 



C. Case #3: Herrmann et al "B-series" 

(g = 1, a = -b, 02 = 62 = 0) 

We next consider the "B-series" simulations of Her- 
rmann et al. Q. This configuration consists of equal- 
mass black holes with equal-magnitude, oppositely di- 



Briigmann et al. noted a sin(2(/> -|- phase) dependenee in S3, 
but concluded that it might bo a non-physical systematic error, 
since they could not find a consistent sin((^ + phase) contribution 
(see their Fig. 13, and corresponding discussion). Our formalism 
explains why the sin(</) -|- phase) term is forbidden, and suggests 
that the oscillation which they observe is probably a real physical 
effect, not a systematic error. 
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rected spins lying in the (e^^^e^^') plane: 



\"3, 

■ (41) 

Though this is obviously a highly symmetric configura- 
tion, it is not one of the special configurations derived in 
Section IIIII and therefore has non- vanishing values for all 
of the final observables / S {m, kj_, A:3, s^, Sg}. As in Case 
^2, the initial spin configuration (|4ip is parameterized 
by two numbers (a and 4>) instead of six {a.i, bi}, so that 
the original spin expansions (|25|) become highly degen- 
erate. Therefore, just as in Case #2, we should remove 
these degeneracies by collecting terms with the same de- 
pendence on a and 0, and rewriting the expansion in the 
form ([57a|) . or equivalently ((37bl [55)) . For convenience. 
Appendix IB 31 provides explicit expressions for the "new" 
coefficients ^(^'J)} in terms of the original spin- 

expansion coefficients J'"i™2ni3|nin2n3^ third order 

in a. 

Since the B-series configuration (|4ip has different sym- 
metries from the superkick configuration (|36|1 . there are 
different rules governing which terms appear in the ex- 
pansions (|37p for each observable. We can discover these 
rules as follows. Applying a parity transformation P to 
the initial configuration (|¥T|) is equivalent to sending (p 
—(p, and hence {sin(j(/)), cos{j(j))} — > {— sin(j(/)), cos(j0)}. 
From this we infer that 

• In the B-series configuration |^ onZy cosine 
(sine) terms can appear in the expansions |5'7[ ) /or 
a scalar (pseudoscalar) quantity f. 

Applying a combined parity and exchange transforma- 
tion PX to the initial configuration (|41|) is equiva- 
lent to the transformation (p (jj + tt, and therefore 
{sm{j(j)),cos{j(j))} {{-ly sin(j0), (-1)^ cos(j</))}. We 
thus infer that 

• In the B-series configuration only terms with 
even (odd) i and j can appear in the expansions 
^37^ for a quantity f that is even (odd) under PX . 

From these simple considerations, a host of interesting 
predictions follow. ^'^ To illustrate this point, consider the 



Herrmann et al. proposed a general kick formula [their 

Eq. (5)] inspired by the Kidder formula for the emission of linear 
momentum [49ll . This formula is linear in the initial spins and 
agrees with ours to this order. A key difference is that they claim 
their formula is only valid when the initial spin components are 
determined at entrance, when the binary reaches the "last" orbit 
or plunge. We claim that symmetry and the inclusion of higher- 
order terms allows our formula to be a valid approximation at 
any initial separation. 



expansions for / G {k^, s^, m}: 
k, = 



'=k|''''a3 + O(a5)]cos(0) 



^kf '^^a=^ + "ki^'^^a^ + Oia')] cos(30) 



'k^^-^^a^ 'k^^''^^a^ + 0{a^)] sin((/)) 
0^ + 0(0^)] sin(3(/)) 



,(3,3)„3 , si,(5,3)„5 



(2'2)„2 , (4,2)^4 



a^ + 0{a^)] sin(20) 



So (4.4) 4 , s„ (6,4) 6 



a'' + 0(a**)] sin(40) 



a 



a-' + O {a")] cos{04>) 



■c. (2,2)^2 , c„ (4,2)^4 



a'*-hO(a^)]cos(20) 

m(2^2)„2^c^(4,2)^4^(5(^6)]^Qg(2(/.) 



(42a) 



(42b) 



(42c) 



(42d) 



(42e) 



Each of these equations makes specific predictions (at 
leading order in a, at next-to-leading order, and so on) for 
how the final-state quantities should vary as a function 
of a and cj). We hope that in the future many of these 
predictions will be tested in detail. 

Currently, there are 7 simulations to which we can 
compare our predictions: the 6 "B-Series" simulations in 
, plus one additional {(p = 0) simulation from an earlier 
work by the same group [T^l . We summarize the relevant 
simulations in Appendix lA 31 Although the authors re- 
port results for the final kick, as well as the final radiated 
energy and angular momentum from each simulation, the 
radiated energy and angular momentum values are not 
accurate enough to constrain the spin expansion beyond 
the trivial zeroth-order (constant) term. Therefore, we 
only consider the final kick. Since the initial spin mag- 
nitude is fixed (a = 0.6 in all 7 simulations) and only (p 
varies, we should rewrite the expansions (|42p in the form 
(j37bp . In particular, the expansions for the final kick 
become 



k^ = "k|^^ cos{(P) + ''k'^' cos{3(p) + 0{a'') (43a) 
fcg =''k^^^ sin((/)) + "4^^ sin(30) + O(a5) (43b) 

As in Case #1, only the magnitudes of the final kicks |k| 
have been published. ^''^ Combining Eqs. (|43ap and (j43b[) . 
we find that |kp = k^ + fc| has the expansion: 

|k|2 = Ao + A2 cos(2(/.) + Ai cos(40) + 0{a^), (44) 



,(3) 



Herrmann et al. Q show the individual components of k in their 
Fig. 10, but without sufficient numerical precision to allow us to 
constrain higher-order terms. 
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Order 


^0 


Ai 




up to c? 


5.04 X 10^ 


-4.23 X 10^ 




up to a'* 


5.05 X 10^ 


-4.20 X 10^ 


-3.92 X 10^ 



TABLE VI: Fits of the kick magnitudes |k|""™ listed in Ta- 
ble |XV] to the fitting formula of . The first column lists the 
a-dependence of the highest-order term appearing in the fit; 
remaining columns provide the best-fit values for the ampli- 
tudes Ai in units of (km/s)^. The first row describes a fit 
to the two second-order terms in Ea. (|44|l . while the second 
row provides best-fit values to all three terms that appear to 
fourth order in a. 



where the amplitudes Ai are given by 



(45a) 



-^^.^l^^vfVk'i^^kf (45b) 



1000 



800 - 



600 



400 



200 




(45c) 



0.6 0.9 1.2 

^(radians) 



1.8 



Note that the expansion only contains terms of the 
form cos(2n0), just as we would expect for a quantity 
with P = -f-l and PX — -|-1. like to or Sg. However, 
since |kp is a composite quantity constructed from the 
more "fundamental" quantities kj^ and , the a° term is 
missing from the coefficient and instead the leading 
order term in A^ is proportional to . This is a crucial 
physical difference between the expansion for |kp and 
the expansions for to and Sg which both contain a non- 
vanishing term. 

When we fit the 7 simulations to the second-order and 
fourth-order forms of Eq. (gll), the x^/d.o.f. is 8.36x 10"** 
and 3.30 x 10"'*, respectively. The corresponding best-fit 
values for the amplitudes Ai are listed in Table IVll and 
the second-order fit itself is displayed in Fig. [5l*^ We see 
that the and A-2, terms already provide an excellent 
fit to the simulations of 0], and there is only marginal 
evidence for the A/^ term. In order to definitively detect 
the presence of the A^ term, additional simulations would 
be needed. 



D. Case #4: Herrmann et al "S-series" 

In this section we consider the "S-Series" simulations 
of Herrmann ei aZ 0] , in which the black holes have equal 



Herrmann et al. Q found that the angle (p (their 9) remains un- 
changed during the final few orbits of inspiral, allowing them to 
apply their kiek formula using the initial spin components rather 
than those at entrance. Their quantities {K^^xi Knaxi Kn^j^} 
correspond directly to our {'^fej,'^ fcj,'' fc^}, leading to values 
Ao = 5.02 x 10^ (km/s)2, Aq = -4.20 x 10^ (km/s)^ consis- 
tent with our results. 



FIG. 5: The kick magnitudes |k| as a function of polar angle 
for the Case #3 configuration given by Eq. (|4H) . The square 
points correspond to the simulations listed in Table IXVl while 
the blue curve shows the second-order fit to Eq. (|44|) with 
amplitudes listed in Table fVll 



mass (q — 1), and spins initially given by 





(46) 



This is the least symmetric spin configuration that we 
have considered so far. Since the configuration is param- 
eterized by two numbers (a and 0) instead of six (a^ and 
fej), we should again use the expansions (j37p . In partic- 
ular, since all 22 simulations in this series have the same 
spin magnitude a — 0.6, we should write the expansion in 
the form (|37bp . Due to the lack of symmetry, most of the 
expansion coeflicients {j'-^K'f^^^} are non- vanishing. In 
Appendix IB 41 we provide relations between these "new" 
expansion coefficients and the original spin expansion co- 
efficients /'"i™2m3|«in2«3 

We first consider the final spin Jl^^i m the e'^' direc- 
tion, and the total energy i?rad emitted in gravitational 
radiation during the inspiral. Both of these quantities 
have expansions of the same form as that of the final 
mass TO, so it is convenient to analyze them in paral- 
lel. As indicated by the fits listed in Table IVIIl there 
is little evidence for terms beyond linear order in a in 
either expansion. These linear-order fits are shown in 
Fig. m These same linear terms, of the form X^^^^^^^b^, 
appeared in our analysis of (anti-) aligned configurations 
(Case #1). They can be understood physically by the 
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Order 



xVd.o.f. 



xVd.o.f. 



^(g.-ad/A/)" 



= (i5.-ad/M) = 



=(£;.-ad/Af)-^ 



0th 



0.321 



0.610 



1.476 



3.45x10" 



1st 



3.96x10" 



0.617 



6.91x10" 



0.0508 



3.66x10" 



8.71x10" 



2nd 



1.29x10" 



0.618 



6.84x10" 



^5.95x10" 



-5.15x10" 



0.0174 



3.68x10" 



8.82x10" 



-1.36x10" 



2.57x10" 



3rd 



9.95x10" 



0.618 



6.85x10" 



-5.78x10" 



iame arguments presented in the final paragraph of that 
Subsection. Though estimates of -Eiad are not available 
or Case #1, we can attempt a quantitative comparison 
jetween the final spins >/|nai/^^^ here and values of S3 
jrovided for that configuration. Equating these expan- 
sions term-by-term at zeroth and first order in a, wc find 



-5.26x10" 



7.31x10" 



-2.34x10" 



0.0161 



SxlO" 



8.82x10" 



-1.41x10 



3.25x10" 



-1.37x10 



3.43x10" 



TABLE VII: Fits for the ^-component of the final black hole 
spin Jina.i/M'^ and radiated energy _Erad/Min Case #4. The 
first column gives the observable being fitted. The second col- 
umn listes the x^/d.o.f. for the best fit. Remaining columns 
provide the best-fit values for the listed parameters. 



0.7 




3 4 

0(radians) 



FIG. 6: The 2— compnent of the final spin Jgnai radiated 
energy -Brad as functions of the polar angle (p between the 
spin b of black hole B and the orbital angular momentum in 
the configuration of Case #4. The square points correspond 
to the simulations listed in Table IXVll while the blue curves 
show the linear fit to Eq. (|37b|l with the first-order coefficients 
listed in Table IVlIl 



^OOOIOOO 1-^0001000-12 

+2s 



(47a) 
(47b) 



" inserting the appropriate values from Tables IIVI and IVIII 
nto the right and left-hand sides of Eq. ((47|) rcspcc- 
— ;ively yields close agreement between '^iJl^ai/^^^)^ — 
'-fc.617 and s0"0|000(to000|000)2 ^ o g26 and between 
(JLJMY = 6.91 X 10-2 and [s™i|ooO(^ooo|ooo)2 ^ 

2s°°°l™°mOOOl™Om™i|ooo]a = 6.96 x lO'^. This agree- 
ment is well within the systematic errors attributed to 
each series of simulations, and suggests that our approach 
of decoupling spin-dependent effects term-by-term shows 
promise. Further simulations will be necessary to deter- 
mine how well this promise is fulfilled. 

We now turn our attention to the black-hole kicks. As 
in the case of the "B-series" considered in the previous 
subsection, our analysis is hampered by only having ac- 
cess to the magnitude of the kicks rather than their in- 
dividual components. In order to most clearly illuminate 
the degeneracies that remain between terms in our gen- 
eral expansion (pS]) . we must unfortunately resort to yet 
another new expansion for this configuration. We expand 
the individual components and squared magnitude as 



k^=;^ cos"(0)rkf ) + ''kf 'sin(0)] , 

00 

fc3=^ cos'"(0)rfcf ' + ^4™)sin(0)] , 

00 

cos"(6')['=A"(")-|-''A'(")sin(6i)] 



(48a) 
(48b) 
(48c) 



m=0 



As previously, we relate the coefficients in this new 
expansion to those in the general expansion in Ap- 
pendix IB 41 This new expansion is useful because to 
third order, '^/cf ^ = "/cf ^ implying that "K^"^ = "^(0) 
and '^^(i) = There are therefore only 2 indepen- 

dent terms in the expansion of |kp when kj^ and ^3 are 
expanded to first order in a. The number of independent 
terms in the expansion of |kp increases to 6 when kj^ 
and are expanded to second order in a, and increases 
again to 10 when the components are expanded to third 
order. The independent coefficients at each order and 
their best-fit numerical values are listed in Table IVIIIl 
while the fits themselves are displayed in Fig. [71 The 
error bars in this case are proportional to the kick mag- 
nitudes |k| themselves unlike in Cases #1 and #2. This 
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1st order 


2nd order 


3rd order 


xVd.o.f. 


19.5 


12.8 


1.50 




6.33 X 10"* 


3.19 X 10"* 


4.87 X 10^ 




— 


-6.16 X lO-* 


-3.03 X 10^^ 




2.75 X 10'' 


7.08 X 10" 


2.92 X 10® 




— 


2.81 X 10^ 


2.53 X 10® 




— 


1.99 X 10^ 


9.99 X 10® 




— 


— 


8.27 X 10® 






-5.36 X lO'' 


-2.95 X 10® 


^7^(4) 






-8.70 X 10^ 








-6.53 X 10® 








-3.05 X 10® 



TABLE VIII; Fits of the squared kick magnitudes |k|^ listed 
in Table IXVll to the fitting formula of Eq.(|lHl). The first 
column shows the coefficients being fitted, while the second, 
third, and fourth columns list the numerical values for these 
coefficients when fits are performed to first, second, and third 
order in a for the individual components of the kicks. The 
magnitudes of the coefficients are given in units of (km/s)^. 




implies that our best-fit spin expansions will naturally 
agree more closely with the smaller values of |k|. 

The angular dependence of the numerically determined 
kicks in this configuration is highly non-trivial, and seems 
to challenge the linear "Kidder" kick formulae contem- 
plated previously in the literature Q . These formulae 
can be reconciled with the numerical results by recogniz- 
ing that they are linear not in the initial spins, but in the 
spins evaluated at merger. Both Kecoii in Eq- (1) of [1] 
and V in Eq. (5) of Q depend on the angles 



=cos-i(£-m«) 



between 



M 



Mb Ma 



(49) 



(50) 



and the components m'*) of an orthonormal triad defined 
at merger. Our triad {e^'^} was defined at the beginning 
of the simulation. For configurations like Case #4 that 
lack high degrees of symmetry, the initial spins appearing 
in S will precess during the inspiral and the orientation 
of the triad {m'^*-'} will vary with respect to a fixed frame. 
This implies that the angles 0^'' appearing in the Kid- 
der formulae implicitly depend on the initial spins, intro- 
ducing non-linearity into the formulae. We attempt to 
capture this non-linear spin dependence explicitly in our 
spin expansions, allowing us to construct fitting formulae 
that only depend on genuine initial conditions. 

The price we pay for making this spin dependence ex- 
plicit is more complicated nonlinear fitting formulae, as 
well as a more cumbersome explicit procedure for relat- 
ing expansions calibrated at different initial stages of the 
inspiral. This procedure for relating different expansions 



FIG. 7: The magnitude |k| of the recoil velocity in km/s as a 
function of the polar angle (f) between the spin b of black hole 
B and the orbital angular momentum in the configuration of 
Case #4. The square points correspond to the simulations 
listed in Table IXVll while the green (dotted), red (dashed), 
and blue (solid) curves show the best fits of Eq. (|48c|) to these 
simulations using terms derived from expanding Eqs. (|48a|) 
and (|48b[) to first, second, and, third order in a respectively. 
The error bars correspond to la errors of 15% in the kick 
magnitudes |k| reported by Herrmann et al @j as the accuracy 
of their simulations. 



will be provided in the second paper of this series. An 
example of how our nonlinear fitting formulae might arise 
from the linear Kidder formulae can be seen in the con- 
figurations of Case #4. According to Eq. (5) of Q, 



fc3 0cS-(i^„m(i'+Xfcm^"J) 



,(2)1 



(51) 



Post-Newtonian expansions [49| reveal that to linear or- 
der the triad {m(*''(a, b)} for equal-mass spinning BBHs 
is related to that in the non-spinning case by 



m 
m 



(1) 

(2) 
,(3) 



(a,b)\ 
(a,b) 
a,b); 




(0,0) 



m(2)(0,0) 
m(3)(0,0) 



(52) 

where 4* is linear in 03 + 63. If we Taylor expand ^' in 
Eq. dSSl) then insert both Eqs. ^ and ^ into the 
Kidder formula ((5T|) . terms prooprtional to 01(03 -I- 63)" 
emerge. We thus see one way in which a formula linear 
in the spins at merger can become nonlinear in the initial 
spins. 

This indeed may be part of the story explaining the 
numerical results in Case #4. The most surprising fea- 



ture of the "S-Serics" kicks is the small kick velocity of 
230 km/s at (/) = 75°. This spin orientation is quite close 
to the "superkick" configuration {(f> = 90°) at which one 
would naively expect the kicks to be maximized. While 
the amplitude of the Kidder formula for is indeed max- 
imized at (/) = 90°, Eq. ([STj) predicts that the superkick 
should have a sinusoidal dependence on the azimuthal an- 
gle at merger as seen in Case #2. This angle will depend 
on the total phase accumulated between the beginning of 
the simulation and merger. The numerical results for the 
S-Series listed in the final column of Table IXVII indicate 
that the duration of the inspiral varies as cos (linear 
in 03 -I- 63). The orbital frequency of equal-mass non- 
spinning BBHs at merger is about uj ~ 0.15M~^ W^ - 
suggesting a final orbital period r = 27r/ijj ~ 42M. This 
estimate of the period is quite comparable to the differ- 
ence in merger times 177.3Af — 138.6A/ — 38.7il/ between 
the successive peaks in |k| a.t (f> = 45° and (f> = 105°. If 
this effect is indeed responsible for the observed kicks in 
Case #4, it helps to explain which third order terms are 
more significant than in previous cases. Further simu- 
lations are necessary to determine if this explanation is 
correct, and how best to account for it in our formalism. 



E. Case #5: Generically oriented initial spins 

Finally, we consider the set of 8 equal-mass simulations 
published in Tichy and Marronetti [ll[. Each simula- 
tion has a different initial spin configuration, and most 
of these configurations have no particular symmetry. The 
initial configurations are not chosen according to a pat- 
tern, but are instead intended to be "generic." In con- 
trast to previous subsections, there is no natural way to 
parameterize them in terms of one or two numbers. This 
means that, again in contrast to previous subsections, we 
cannot use symmetries and degeneracies to significantly 
reduce the "effective" number of coefficients in the spin 
expansion. We must therefore truncate our spin expan- 
sions at an order for which the number of independent 
terms is fewer than the number of simulated data points 
if we hope to non-trivially test our formalism. 

For each of the 8 simulations in [ll[ , the authors quote 
the final kick magnitude |k|, the final spin magnitude |s|, 
and the final mass m. Is this enough information to test 
the spin expansion formalism? The answer is "no" for 
|k|, and "yes" for |s| and m. In the previous subsec- 
tions, we have seen that in the general case we should go 
to second- or even third-order in the spin expansion to 
achieve a good fit for |k|; but already at second-order, the 
independent coefficients in the general spin ex pan sion for 
|k| outnumber the 8 values of |k| provided by [il|. Thus, 
we cannot use the |k| data to perform a non-trivial test. 
The situation for m (and even for |s|) is better because, 
as we have seen in previous subsections, the linear terms 
in the spin expansion seem sufficient to provide a good 
fit to the data (except in special symmetric cases where 
these linear terms vanish, so that the second-order terms 
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Quantity 


Results 


m 

Eq. If53a)l 


xVd-o.f. = 0.275 
free parameters 


1^1 

Eq. JSl 


xVd.o.f. = 0.671 

free parameters 


is| 

Eq. l(53d)) 


xVd.o.f. = 0.0707 
3 free parameters: 
|si°°l''°''| « 0.1871 
|sr""V 0.2086 

cose ^ -0.974 



TABLE IX: Fits to the final masses m and spin magnitudes 
|s| of the simulations listed in Appendix lA 51 The first row 
lists the x^/d.o.f. of fits of Eq. (|53a|) : we use the values of 
the coefficients determined in Case #1 leaving no free pa- 
rameters. The second row shows fits of Eq. (|54|) to |s|; using 
the appropriate coefficients determined from Case #1 again 
leaves no free parameters. The third row fits the full formula 
of Eq. H53d|l to jsj; we list the x^/d.o.f. and best-fit values of 
the 3 new parameters. 

become important). Let us test whether this simple be- 
havior continues to hold for the generic initial spin con- 
figurations of (llj . At linear order, the mass and spin 
magnitude are given by: 

m =m™°l"°" + mO"i|°°"(a3 + 63), (53a) 

100|000/ , I \ , 010|000/ , , N /rou^ 

s_L = Sj^ (ai + bi) + Sj_ (024-02), (53b) 

OOOlOOO , 0011000 / , 7 ^ /ro N 

53=53 (03 + 03), (53c) 

M =^K\' + sl (53d) 

The four coefficients {m°°°\°°°, ^ooi|ooo^ gOOO|ooo^ ^ooi|ooo^ 
were already accurately determined by the large set 
of simulations in Case #1. We therefore fix these 
coefficients to the values listed in Table IIVI Eq. (|53a[) 
thus predicts the final mass m with zero free parameters. 
If we Taylor expand the square root in Eq. (j53dp for |s|, 
keeping terms up to linear order in the initial spins, we 
obtain 

I I OOOlOOO , 0011000/ , u \ ^^;/l^ 

Nl-Sg +S3 ' (03 + ^3) (54) 

which again has no free parameters. Alternatively, if we 
do not Taylor expand the square root in Eq. (|53dp , then 
the expression for |s| has three free parameters: the mag- 

, r .1 • J- loolooo J oiolooo , 

nitudes of the coefhcients Sj_ ' and s_|_ ' and the 

angle 8 between them. In Table HXl we list the x^/d.o.f. 
and the best-fit values of any free parameters from fitting 
to Eqs. ((55a|) . ([5i)) . and (|53d|) . The fits themselves are 
shown in Fig. [S] 

The goodness of these fits speaks for itself. The most 
remarkable result is that Eqs. (|53ap and ((54|) — which 
have no free parameters and only depend on 03-1-63 — 
do an excellent job (red points in Fig.[S]). For all their im- 
portance in generating the "superkicks," the components 
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0.98 



0.97 - 




0.6 - 



0123456789 

simulation number 

FIG. 8: Final masses m and spin magnitudes |s| for the 8 
simulations Ai of Case #5 listed by simulation number i. The 
black points show the numerically determined values and la 
error bars listed in Appendix I A 51 The red points in the top 
and bottom panels show the predictions of Eqs. (|53a|l and 
(|54l) respectively. These use the best-fit values of coefficients 
obtained in Case #1 and fit no additional parameters. The 
blue points in the bottom panel show the best-fit of Eq. (|53d[l 
to |s|, with coefficients listed in Table HXl 

of the initial spins in the orbital plane {a^^ bi, 62} seem 
to have little effect on the final masses and spins. Only 
one simulation (#5) is fit poorly by the red points; not 
surprisingly it is the configuration with the largest spin 
projection in the orbital plane. The red points over- 
predict m by failing to account for the energy carried 
away by the gravitational radiation sourced by the pla- 
nar spins, and underpredict |s| by neglecting the contri- 
bution of sx. These results are provocative, but remain 
provisional until verified by further simulations. 



V. CALIBRATING THE SPIN EXPANSION 
WITH NEW SIMULATIONS 

In this section, we suggest a relatively small set of simu- 
lations (10 equal-mass simulations, and 16 unequal-mass 
simulations) with initial spin configurations specially cho- 
sen to allow all of the spin-expansion coefficients up to 
second order to be determined uniquely. Once the spin 
expansion is calibrated in this way, it becomes fully pre- 
dictive — i.e. it predicts the simulated observables (to 
second-order accuracy), for any initial spin configuration. 

This section is organized into four subsections. Sub- 



section explains why currently available simulations 
are not sufficient to calibrate the spin expansion. Sub- 
sections IV Bl and IV CI suggest an explicit choice of 10 
equal-mass and 16 unequal-mass simulations with ini- 
tial spin configurations suitable for this calibration, and 
provide corresponding formulae for the spin expansion 
coefficients in terms of the results of these simulations. 
Subsection IV Dl collects several additional comments for 
readers who are interested in pursuing or extending the 
program suggested here. The importance of performing 
these simulations is discussed later in Sec. I VII 

In Appendix [Dl we generalize the calibration proce- 
dure to take account of the inevitable systematic errors 
arising from the simulations themselves and from trun- 
cating the spin expansion at finite order. In particular, 
given N simulations of known covariancc, we derive the 
minimum-variance unbiased estimators for the spin ex- 
pansion coefficients, as well as the corresponding covari- 
ance matrix for these estimators. 



A. Breaking degeneracies 

In the previous section, we examined in detail the pre- 
dictions of the spin expansion for the 5 different con- 
figurations (Case #1 through Case #5) that have been 
simulated to date. Although at first glance it might seem 
that the total number of currently available simulations 
is more than sufficient to calibrate all of the expansion 
coefficients up to second order, in practice we had to 
perform new fits for each of the 5 configurations. These 
new fits were required for three distinct reasons. Firstly, 
many published works have only provided the final kick 
and spin magnitudes |k| and |s|, not the individual com- 
ponents fcj and Sj. Combining the spin expansions for 
individual components to predict final magnitudes intro- 
duces degeneracies between terms. 

Secondly, many of the available simulation sets study 
the same highly symmetric configurations — like spins 
aligned with the orbital angular momentum (Case #1), 
or the superkick configuration (Case #2). These par- 
ticular configurations have justifiably garnered much of 
the attention, both because of their astrophysical inter- 
est and because they cleanly illustrate several key fea- 
tures of spinning BBH merger. Nevertheless, since the 
initial spin components in these configurations are either 
purely aligned with the orbital angular momentum (Case 
#1) or purely perpendicular (Case #2), they fail to con- 
strain the coupling between aligned and perpendicular 
spin components allowed by symmetry beyond linear or- 
der. We found that such terms were required to attain 
a good fit to the kicks of Case #4, suggesting that fur- 
ther simulations are indeed necessary to determine the 
coefficients of these terms. 

Thirdly, different groups performed their simulations 
at different values of the initial dimcnsionlcss orbital sep- 
aration r/{Mg^ + M^). Recall however from Sec. [Ill that 
our spin-expansion coefficients are defined at a fixed value 
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w 


X 


y 


z 


p 


+ 1 


+1 


-1 


-1 


X 


+ 1 


-1 


+1 


-1 


e.g. 


m, S3 


ki, k2 


Sl, S2 





TABLE X: Eigenvalues under P and X for the generic final 
quantities w, x, y, and z, in subsection IV Bl along with ex- 
amples of physical quantities with these transformation prop- 
erties. 



of r/{M^ + Mb) — or, more generally, at a fixed value 
of the dimensionless inspiral parameter ip. While post- 
Newtonian techniques may be able to relate the spin- 
expansion coefficients at different values of ^p, doing so 
will add an additional layer of complication and poten- 
tial systematic error to the calibration process. We treat 
this issue in a forthcoming paper. 

As existing simulations are insufficient to fully con- 
strain the spin expansion, in the following subsections 
we explicitly provide the initial spin configurations for a 
small number of new simulations with which we will be 
able to achieve the desired calibration. Since our spin- 
expansion coefficients remain functions of the mass ratio 
5, a new set of simulations is required in principle at each 
value of q. In practice, we hope that the g-dependence 
of our coefficients is sufficiently smooth that we can in- 
terpolate between coefficients calibrated at a small set of 
mass ratios. 



B. Equal-mass (g — 1) BBH mergers 



Final quantities for equal-mass {q = 1) BBH mergers 
are characterized by their eigenvalues ±1 under parity 
P and exchange X. In this subsection, we will use the 
4 variables {w,x,y, z} to denote generic final quantities 
with the 4 possible combinations of these eigenvalues as 
shown in Table 1x1 

To second order in the spin expansion, w, x, y, and z 



eqO: 


none 


eql; 


eq2: 02 


eq3; ai, 


eq4: 61, 03 


eq5; a^, Oj, 


eq6: aj, 62, ^3 


eq7: ai, &i 


eq8: aj, &2 


eq9: a^, 


N.B. afV4^) 



TABLE XI: A suggested set of 10 simulations (denoted eqO 
through eq9) which can simultaneously determine all of the 
coefficients up to second order in the general spin expansions, 
Eqs. (|55p . for equal-mass {q = 1) BBH merger. For each 
simulation, we list the non-vanishing initial spin components. 
These spin components may all be chosen independently apart 
from the requirement that the values of as differ between sim- 
ulations eq3 and eq4 {a'^^ 7^ ^s*'')- While many different sets 
of configurations satisfy this requirement, it is particularly 
convenient to choose the set given by Eq. (|56|l . 



are given by: 

+ i«ooi|o"0(a3 + &3) + u;002|oo"(a2 + 62) 
-f U.200I000 (^2 + 62) + ^0201000 ^al + bl) 

+ i;;""|0"0(aia2 + ¥2)+^'°'"°'°(«i^2+^a2) 
+ wi°"|i"°ai6i+u;"i"l"i°a2&2-Hu;°°i|™ia353 (55a) 
a; = x0°i|™°(a3-53)+x"02|000(„2_^2) 

^^200|000(„2_^2)^,^020|000(^2_^2) 

+ x""l™°(aia2-&i62)+xi™l°i°(ai52-5ia2) (55b) 
y = yi°°l°°°(ai+6i)+2/"i°l°°"(a2 + 52) 

+ y'°"\""'{a,b, + b,a,) + y'>''>\''\a,b, + b,a,) (55c) 

z = zioo|ooO(a^-6^) + ;^oi"|0"0(a2-62) 
-fzi°i|°°°(a,a3-6i63) + 2""l°"°(a2a3-5263) 
+ z''°\'°\a,b,-b,a,) + z"'°\''^\a,b,^b,a,) (55d) 

Note that the expansions for x, y, and z each contain 6 
coefficients, while the w expansion has 10 coefficients. 

In Table |Xll we suggest a set of 10 simulations de- 
signed to determine all the coefficients of Eq. (|55p at the 
minimum computational cost. We name these 10 configu- 
rations {eqO, . . . , eq9}. The rth simulation (r = 0, . . . , 9) 
has initial spin components {a^\ ^i''}; ^i^d leads to final- 
state values w = w^'"-', x = x'--'^\ y = y^'^\ and z = 

In each of the proposed simulations, only those initial 
spin components listed in Table IXII are non-zero. Al- 
though all 10 of the simulations {eqO, . . . , eq9} are neces- 
sary to determine the w expansion coefficients, only the 6 
simulations {eql, . . . , eq6} are required to determine the 
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X, y, and z coefficients. Thus, if one judges that a 10- 
simulation set is too expensive for one's computational 
budget, the 6 simulations {eql, . . . , eq6} provide a less 
ambitious but still useful alternative (for example, they 
fully calibrate the expansion for the kick components k^). 

The simulations in Table IXll lift all degeneracies up to 
second order in the initial spins. Any choice for the values 
of the non- vanishing initial spin components {a-\^\b[^''} 
in this table will yield a unique solution for all of the coef- 
ficients, provided Og^'' ^ ^s^''- We can use this additional 
freedom to make the final expressions (|57p for the ex- 
pansion coefficients in terms of the simulated observables 
{'w^''\ x^'^\ y^^\ z^^''} as algebraically simple as possible. 
In particular, if we choose the initial non-vanishing spin 
components in Table IXj to satisfy 



"3 
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(56) 



then Eqs. (|55p can be inverted to yield the simple result: 
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/a2'^3 

^(6)„^(3)+^(2)]/^^Q^^ 



where we have defined w'^^'^ = w^''^ — y;000|ooo. 

Once Eqs. ([55|) have been calibrated with the 10 simu- 
lations {eqO, . . . , eq9}, they can predict the results of any 
additional equal-mass simulations to second-order accu- 
racy. 



C. Unequal-mass (g 7^ 1) BBH mergers 

Calibrating the spin-expansion coefficients for unequal- 
mass (g 7^ 1) BBH mergers proceeds similarly to the 
equal-mass case discussed in the previous subsection. 

Without loss of generality, we can choose < g < 1, 
since exchange symmetry X relates coefficients for mass 
ratios q and l/q. However, at a fixed value of q (7^ 1), 
exchange X no longer restricts the terms appearing in 
the expansions. So, for the purposes of this section, 
there are only two types of final quantities: scalars 
and pseudoscalars. Scalars (like {m, fcj^, fcj, S3}) will be 
represented by the variable u, while pseudoscalars (like 
{s]^, S2, A:3}) will be represented by v. 

To second order in the spin expansion, u and v are 
given by 



000|000 



000 020 h 2 



^^200|000^2^^000|200^2^^020|000„2^^i 
+ W""l|"""a3+7.°"°l°°l63+7.°"2|000^2^^^000|002^2 

+ 7."0|"00aia2+«°''°'""&i^2+^'"°'"'"ai&2+"°'°"°°&ia2 
+ ^zi°°|i°°aA+u"i"l°i%62+ii"°i|""ia363, (58a) 



lOOfOOO 



OOOlOlOi 



v = v ai+i>°"°|i°"&i+«"i"l°"%+«' 

+ t;ioi|o"Oaia3+i;0"o|ioi6i63+i;0"loo"a2a3-Fi;0"ol""6263 
+ i;i°°l™iai63+«°°i|^°°&ia3+z;°i°l°°ia263+z;OOi|"io&2a3. 

(58b) 

Note that the expansion for u contains 16 coefficients, 
while the expansion for v contains 12 coefficients. 

In Table IXIIl we suggest a minimal set of 16 simu- 
lations needed to simulataneously determine all of the 
coefficients in Eqs. ([55)1 . We have named these con- 
figurations {uneqO, . . . , uneql5}. As in the equal-mass 
case, the rth simulation (r = 0, . . . , 15) has initial spin 

components {af\b'f^}, and leads to final simulated ob- 
servables u = u^^^ and v = v^'^\ In each of the 
proposed simulations, only the initial spin components 
listed in Table IXIIl are non-zero. While all 16 of the 
simulations {uneqO, . . . , uneql5} are required to deter- 
mine the u expansion coefficients, the 12 simulations 
{uneql, . . . , uneql2} suffice to calibrate the v coefficients. 
If one is only interested in final quantities with odd par- 
ity (like si, S2, and fcs), one can reduce the computing 
time by ~ 25% by only performing the 12 simulations 
{uneql, . . . , uneql2}. 

The simulations in Tablc lXlIl vicld a unique solution for 
all of the coefficients in Eqs. as long as one chooses 
flg^-* 7^ af^ and Jj'p 7^ . As before, one can choose 
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uneqO: none 


uneql: 


uneq2 


"2 


uneq3: 61 


uneq4 


62 


uneq5: a^, 03 


uneq6 


"2, as 


uneq7: 61, & 3 


uneqS 


62, 63 


uneq9: a^, bi, 03 


uneqlO: 02, ''2, "3 


uneqll: Oi, 62, &3 


uneql2: 61, aj, 63 


uneql 3 


fli, 02 


uneql4 


62 


uncql5 


13, ^3 


N.B. arVaf,ferV;>f 



TABLE XII: A suggested set of 16 unequal-mass {q 7^ 1) sim- 
ulations that constitute a minimum set necessary to calibrate 
the spin-expansion coefficients in Eq. (|58|l . For each simula- 
tion, we list the non-vanishing initial spin components. These 
spin components may all be chosen independently, apart from 
the constraints a'^^ / 03^' and 63^' 7^ feg*', but it is particu- 
larly convenient to choose components that satisfy Eq. (|59|) . 
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(60b) 



where we have defined u^^'^ = u^''-' — ^0001000. 

Once this calibration has been achieved, Eqs. ([55]) will 
predict the simulated observables to second order accu- 
racy, for any initial spin configuration. 

D. Technical points 



the values {al ,6,- } appearing in this table to make the 
inversion of Eqs. (|58p particularly simple. If these initial 
spin components satisfy 



a. 



'2 



(5) 



(10)_ 



-_;.(8)_ 



/92 = 
/?3^ 



+b, 

+b. 



(4) 
2 " 
(7)_ 



,(13) 
^1 

,(13) 
,(15) 



(59) 



,(14) 



the inverted equations for the spin-expansion coefficients 
take the comparatively simple form 
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000 
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^001 
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-(S(6)_u(2) 
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y002 


000 


^[(y(5)_y(l)) 


+ (u(6)_u(2) 


)]/2ai 


yOOO 
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= [(y(7)_y(3)) 


_(u(8)_s(4) 


)]/2/33 


^000 
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= [(?l(7)-u(3)) 


+ (y(8)_y(4) 


)]/2/3| 


yOOl 


001 


= [y(15)_u(7)„ 


u(5)-Hm(3)+ 
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In this final subsection, we collect a few additional re- 
marks that will be of interest to readers who wish to 
pursue or extend the program suggested in this section. 

So far, we have discussed the initial spin orienta- 
tions necessary for calibrating the expansion coefficients. 
What about their absolute magnitudes |a| and |b|? It 
is best to choose the initial spins to be rather small for 
two reasons. Firstly, the errors introduced by neglect- 
ing terms beyond second order are fractionally smaller 
for small initial spins. The values of the coefficients 
given by Eqs. (|57)) and ([5U)) will therefore be closer to 
their true, unbiased values. The second point relates 
to the way in which initial conditions for simulations 
are presently specified. Most groups currently take the 
3-metric jij on the initial 3-dimensional spatial hyper- 
surface to be conformally fiat. While an isolated, non- 
spinning (Schwarzschild) black hole has conformally flat 
spatial hypersurfaces, a spinning (Kerr) black hole does 
not (ssj and neither do BBHs (spinning or non-spinning). 
As the initial spins increase, choosing the initial 7^ to be 
conformally flat is expected to become an increasingly 
poor description of realistic BBH initial data. This choice 
will therefore lead to correspondingly larger systematic 
errors in determining the coefficients in the spin expan- 
sion. The initial spins should be chosen small enough 
to minimize these problems, yet large enough that the 
second-order effects we are seeking are not swamped by 
other systematic errors in the numerical simulations. 

As these systematic errors are inevitable, it may be 
fruitful to reinterpret the coefficients on the left-hand 
sides of Eqs. ([57)1 and ([50]) . Instead of regarding these co- 
efficients as the algebraic solutions to Eqs. ((55)) and (|58)) . 
we consider them to be estimators (t()™i™2ni3|nin2n3 ^ 
-mimsmalnmsna^ . . . ) constructed from the simulated ob- 
servables {w^'^y x'^'^\ . . .} and the estimated initial spin 
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components {a^ , }■ Once errors are taken into ac- 
count, it may be desirable to use larger sets of simulations 
to beat down the noise associated with our estimators. 
We pursue this approach in Appendix [Dl where we de- 
rive minimum-variance estimators constructed from TV 
simulations, and provide the covariance matrix for these 
estimators in terms of the covariance matrices associated 
with the estimated observablcs {w^"^^ , x*^'"^ , ■ ■ •} and initial 
spin components {a^^^b'f^}. 

In the future, it may be interesting to extend this sec- 
tion's second-order calibration up to third order. Al- 
though all of the third-order contributions to the mass 
TO and spin s identified in Section IIVI were small cor- 
rections, the final kicks k in Case #4 seemed to exhibit 
considerable third-order effects. Recoils for generic ini- 
tial spin orientations have not yet been adequately sim- 
ulated to determine whether these third-order effects re- 
flect genuine physical behavior or are merely artifacts of 
systematic errors within the numerical codes. If spin ex- 
pansions calibrated to second order according to the pro- 
gram outlined in this section fail to describe BBH merg- 
ers with generic initial spin orientations, we may want 
to test whether a third-order expansion can remedy ob- 
served discrepancies. Calibrating to third order will re- 
quire 12 additional simulations (for a total of 104-12=22) 
in the equal-mass {q = 1) case, and 28 additional sim- 
ulations (for a total of 16-f 28=44) in the unequal-mass 
(g ^ 1) case. In Appendix [C] we have provided explicit 
third-order expansions of the 4 variables {w, x,y, z} in 
the equal-mass case. These expansions can be inverted 
to obtain formulae for the third-order coefficients similar 
to the second-order inversions of Eqs. ((55|) . These formu- 
lae can then be used to identify an optimal choice of 22 
simulations from which all coefficients up to third order 
can be calibrated. 



VI. DISCUSSION 

In this paper, we have developed and tested the "spin 
expansion" formalism. This is the following simple idea. 
Long after merger, we regard any final (dimcnsionlcss) 
quantity / — such as the kick velocity k or spin vector s 
of the final Kerr black hole — as a function /(■)/), a j, 6^) 
of the 8 "initial" (dimensionless) quantities {ip,q,a^,b.^} 
necessary to specify the initial configuration of a BBH 
in circular orbit. Then we Taylor expand this function 
around = = 0, and use three symmetries (rotation 
R, parity P, and exchange X) to significantly restrict 
the terms that can appear in the expansion. Finally, we 
interpret the leading-order terms in the Taylor expan- 
sion as leading-order predictions for /, while the next- 
to-leading terms in the expansion arc the ncxt-to-leading 
predictions, and so on. 

To us, it seems genuinely surprising that the final state 
of the complicated non- linear process of binary black hole 
merger can be usefully described by such a simple-minded 
approach. This simplicity should be regarded as another 



discovery which has come from the recent breakthroughs 
in numerical relativity. 

In the Introduction to this paper, we listed some of 
the advantages — both practical and conceptual — of the 
spin expansion formalism. It may be helpful to look back 
at this list, now that we have had a chance to introduce 
and explore the formalism in detail. Here we would just 
like to highlight three new discoveries which came from 
applying the spin expansion to simulations in Sec. IIVI 
and which illustrate the potential of this approach. 



A. Three highlights 

First, we have discovered a new third-order spin de- 
pendence of the kick velocities in the "superkick" con- 
figuration considered in Section llVBI These third-order 
modulations, clearly revealed in Fig. 51 are present in 
the simulations of p, 0] but went unnoticed because with 
amplitudes less than 100 km/s they are dwarfed by the 
primary linear superkicks. We were able to find them 
because the spin expansion made a specific prediction 
for the next-to-leading contribution: it told us to look 
for a contribution to proportional to with triple 
the fundamental (linear) superkick frequency. Empirical 
fitting formulae — linear in spins, and inspired by post- 
Newtonian results — provide acceptable fits for these su- 
perkick simulations, but our discovery shows that there 
is more to learn if one is willing to go beyond these fitting 
formulae. 

Second, we have discovered a new second-order spin 
dependence of the radiation energy ii^rad, the radiated 
angular momentum Jrad and the final spin Sg in the su- 
perkick configuration. The spin expansion predicts that, 
since these three quantities {-Brad, Ji-ad, S3} are all charac- 
terized by the same transformation properties (P = -|-1, 
X = +1), they should all exhibit the same next-to- 
leading-order behavior: A + Bcos{2(j> + phase). Again 
this behavior is present in the simulations of p, 01 , and 
is clearly displayed in Fig. [3l but without the guidance 
of the spin expansion, it went unnoticed in [dl, and was 
dismissed as a possible numerical artifact in [7|. 

Third, we wish to highlight the remarkable agreement 
between the predictions of the spin expansion and the 
simulations of generically oriented spin configurations 
fu\ considered in Sec. IIVE] (Case #5). This agreement is 
illustrated in Fig. [8l where the black points are the sim- 
ulations results and the red points are the predictions. 
We emphasize that the red points are genuine predic- 
tions — i.e. there were no free parameters in these fits, 
since all of the relevant coefficients had already been cali- 
brated by the simulations in Case #1. The red and black 
points only disagree for one of the 8 simulations in Case 
#5 and, as explained in Sec. IIVEI this disagreement is 
easily understood. So far, mumcrical relativists have fo- 
cused mostly on highly symmetric configurations like the 
aligned case in Section IIV Al and the superkick configu- 
ration of Section flV Bl This is probably because of the 



24 



expected complications from non-linear spin precession 
in the generic case. Herrmann et al. Q observe these 
precessions in their "S-Series," and note that they make 
it impossible to use the post-Newtonian-inspired fitting 
formula for the final kicks in this case. Fig. [8] seems to 
provide evidence that our spin expansions continue to 
apply, even in the presence of these precession effects. 



B. Future directions 

Let us end by briefly mentioning a few directions for 
further study. 

First, it would be extremely fruitful to calibrate the 
spin expansion coefficients, up to second or third order. 
As explained in Sec. |Vl currently available simulations 
leave many degeneracies among spin expansion coeffi- 
cients, even at first and second order. To rectify this 
problem, in Sec.|V]we suggest a small set of simulations 
— 10 equal-mass simulations and 16 unequal-mass simu- 
lations — and explicitly show how these would uniquely 
determine all of the spin-expansion coefficients up to sec- 
ond order. Once these coefficients are calibrated in this 
way, the spin expansion becomes fully predictive: given 
any initial spin configuration, it predicts the final results 
{m, /Cj, s,j} with second-order accuracy. In addition to fa- 
cilitating tests of the spin expansion, it is clear that this 
result — a set of simple formulae which predict the final 
state of BBH merger given the initial state — would be 
of enormous interest from the standpoint of astrophysical 
and cosniological applications. For example, our spin ex- 
pansion precisely encapsulates the relevant information 
for incorporating the recent discoveries of numerical rel- 
ativity into cosmological simulations of BBH merger in 
the context of structure formation. It would also be inter- 
esting from a purely theoretical standpoint. The initial 
spin configurations we have identified in Section |V] pro- 
vide a systematic approach for seeking qualitatively new 
behavior in the unexplored regions of BBH parameter 
space. Any unexpected constraints, patterns, or rela- 
tionships among the calibrated coefficients (beyond the 
ones we have used thus far in our construction) could 
indicate interesting new dynamical effects or symmetries 
of the system. 

Second, we have mentioned that a final quantity / 
may be regarded as function on an 8-dimensional space 
{■0, q, Oj, bi}. The spin expansion elucidates the structure 
of the 6-dimensional subspace parameterized by {a^, b^}, 
but we would also like to know the behavior along the 
■0 and q directions. In a follow-up paper, we consider 
how post-Newtonian techniques may be used to to ex- 
plore the T/j-dependence of the spin expansion coefhcients. 
Some insights may be gained by an analogy with effective 
field theory, where the renormalized coupling constants 
depend on the momentum scale at which they are de- 
fined, although the physical predictions of the theory do 
not. Determining the g-dcpendcncc of the spin-expansion 
coefficients |24| seems less straightforward, and is an in- 



teresting topic for future research. 
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APPENDIX A: TABLES OF SIMULATED FINAL 
KICKS, SPINS, AND MASSES 

1. Case #1: g = 1, = = 

We use the 28 simulations of [l^ with non-zero re- 
coils to test our expansions for the final kicks in the 
case of equal- mass (g = 1) binary black holes with spins 
aligned (or anti-aligned) with the orbital angular momen- 
tum. This is the largest and most recent series of simu- 
lations for this configuration published at the time this 
manuscript was prepared. The numerical estimates of 
|kx| are avilable in Table 1 of [lB|. We adopt their la 
errors of 8 km/s for the kick magnitude. 

We constrain our expansion for the final spin S3 by 
performing a joint fit to the full set of 38 simulations in 
[Tst and the 10 simulations of [l3l that they consider to be 
relatively free from the numerical dissipation of angular 
momentum. We use the proposed Icr errors of 0.01 and 
0.02 for the simulations of [l3| and [l3l respectively. 

To constrain our expansion for the final mass m, we 
perform a joint fit to 3 series of simulations [i2,[i3j[I3j as 
each series consists of only a small number of individual 
simulations. These simulations are summarized in Table 
IXIIII For points Ai through A4 we assume fractional 
errors on the radiated energy identical to those provided 
in [l^l for the final kicks. For points Bi through Bg 
we assume errors on the final masses of 0.5% of the ini- 
tial energy Madm, as [l3| only claims to conserve energy 
to this accuracy. Finally, we use the highest-resolution 
simulations of all 11 initial data sets listed in Table I 
of [l3|. While they claim that numerical dissipation of 
angular momentum makes estimates of S3 unreliable for 
031^3 > 0.75, the final masses are largely unaffected as 
shown in their Fig. 4. We therefore use both simulations 
with 03,63 > 0.75, and assume for all simulations errors 
on m of 0.004 consistent with their claimed resolution 
limits. 





as 


&3 




^1 


0.2 


-0.2 


0.9526 ± 0.0023 




0.4 


-0.4 


0.9521 ±0.0017 


^3 


0.6 


-0.6 


0.9519 ±0.0014 


Ai 


0.8 


-0.8 


0.9521 ± 0.0028 


Bi 


0.584 


-0.584 


0.9536 ± 0.0049 


B2 


0.584 


-0.438 


0.9507 ± 0.0049 


Bz 


0.584 


-0.292 


0.9482 ± 0.0049 


Bi 


0.584 


-0.146 


0.9461 ± 0.0049 


Bs 


0.584 





0.9439 ± 0.0049 


Be 


0.584 


0.146 


0.9412 ± 0.0049 


Br 


0.584 


0.292 


0.9376 ± 0.0049 


Bg 


0.584 


0.438 


0.9344 ± 0.0049 


B9 


0.584 


0.584 


0.9315 ± 0.0049 


Ci 


-0.90 


-0.90 


0.970 ± 0.004 


C2 


-0.75 


-0.75 


0.968 ± 0.004 


C3 


-0.50 


-0.50 


0.963 ± 0.004 


C4 


-0.25 


-0.25 


0.958 ± 0.004 


Cs 


0.0 


0.0 


0.951 ±0.004 


Ce 


0.25 


0.25 


0.944 ± 0.004 


C7 


0.50 


0.50 


0.933 ± 0.004 


Cs 


0.62 


0.62 


0.926 ± 0.004 


C9 


0.75 


0.75 


0.916 ±0.004 


Cio 


0.82 


0.82 


0.909 ± 0.004 


Cii 


0.90 


0.90 


0.906 ± 0.004 



TABLE XIII: Final mass data for Case ^1: equal-mass {q — 
1) binary black holes with spins aligned (or anti-aligned) with 
the orbital angular momentum. Points Ai through A4, are 
from [l^ . points Bi through _Bg are from [l^ . and points Ci 
through Cii are from [itI ]. 
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a 




7 num 
^3 


Jrad /M^ 


%-Brad 




0.515 


1.571 


1833 ± 30 


0.248 ± 0.003 


3.63 ±0.01 


A2 


0.515 


0.785 


1093 ± 10 


0.244 ± 0.003 


3.53 ±0.01 


A3 


0.515 


3.142 


352 ± 10 


0.246 ± 0.004 


3.57 ±0.01 


Ai 


0.515 


4.712 


-1834 ±30 


0.249 ± 0.003 


3.63 ±0.01 


As 


0.515 


3.304 


47± 10 


0.245 ± 0.005 


3.55 ±0.02 


Ae 


0.515 


0.0 


-351 ±10 


0.246 ± 0.003 


3.57 ±0.02 




a 




7 num 
^3 


S3 


— 


Bi 


0.723 


0.0 


2680±94 


0.6859±0.0009 




B2 


0.723 


0.524 


2310±94 


0.6856±0.0009 




Bi 


0.723 


1.047 


1150±94 


0.6897±0.0009 




Bi 


0.723 


1.571 


-340±94 


0.6931±0.0009 




Bs 


0.723 


2.094 


-1590±94 


0.6934±0.0009 




Be 


0.723 


2.618 


-2400±94 


0.6892±0.0009 




B7 


0.723 


3.142 


-2690±94 


0.6859±0.0009 




Ba 


0.723 


3.665 


-2320±94 


0.6856±0.0009 




B9 


0.723 


4.189 


-1160±94 


0.6897±0.0009 




Bio 


0.723 


4.712 


320 ±94 


0.6932±0.0009 




fill 


0.723 


5.236 


1580±94 


0.6935±0.0009 




B12 


0.723 


5.760 


2400 ±94 


0.6892±0.0009 





TABLE XIV: Data for Case #2: equal-mass [q = 1) binary 
black holes with equal and opposite spins in the orbital plane. 
Points Ai — Ad are from while points Bi — B12 are from 
O. These papers present final quantities different from those 
discussed in this paper, but whose spin dependence can be 
readily analyzed in our formalism. The scalar Jrad is the to- 
tal angular momentum radiated in gravitational waves, listed 
in units of where M = Ma + Mb is the sum of the horizon 
masses of the initial binary black holes. %i5rad is the percent- 
age of the initial energy radiated in gravitational waves. 



3. Case #3: Herrmann et al "B-series" 



2. Case #2: 5 = 1, ax = -bx, 03 = 63 = 



The data points for this case are summarized in Ta- 
ble IXIVI Campanelli et al. provides error estimates 
for their final observables, which we assume represent 
true 1(7 statistical error bars. Briigmann et al. does not 
provide error estimates for the individual simulated data 
points Bi, however they do claim 95% confidence lim- 
its of ±2% on their maximum kick amplitude of 2, 725 
km/s and ±5 x 10~* on their mean spin ao = 0.6891 
as determined from the black-hole ringdown. We as- 
sume these values correspond to 2a error bars on each 
parameters, and that they were derived from 12 inde- 
pendent data points. This leads to crude Icr errors of 
1/2 X 712 X 0.02 X 2,725km/s = 94km/s on each kick 
and 1/2 X ^/V2 x 0.0005 = 0.0009 on each final spin. 



The data for this case are summarized in Table IXVl 



4. Case #4: Herrmann et al "S-series" 

The data for this case are summarized in Table IXVH 



5. Case #5: The generic case (Tichy-Marronetti) 

The data for this case are summarized in Table IXVIII 
Even at second order in the initial spin magnitude a, 
there are too many independent non-degenerate coeffi- 
cients to fit with only 8 simulations. We therefore only 
attempt to fit the final masses Mf/M and spin magni- 
tudes JfjyP' as these can be fit with linear order terms 
in our formalism. We assume errors on these quantities 
that are 20% of the radiated energy (Af^°^ - Mf)/M 
and radiated angular momentum [J^^ — Jf)/M'^. 
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a 


<!> 


|k|"™ 




-Erad / M 


A 


0.6 


0.0 


285±12 


— 





Bi 


0.6 


0.349 


427 ± 64 


0.24±0.036 


0.033±0.005 


B2 


0.6 


0.524 


544 ± 82 


0.24±0.036 


0.033±0.005 


B-i 


0.6 


0.873 


761 ± 114 


0.25±0.038 


0.034±0.005 


Ba 


0.6 


1.222 


908 ± 136 


0.25±0.038 


0.034±0.005 


B5 


0.6 


1.396 


945 ± 142 


0.25±0.038 


0.034±0.005 


Be 


0.6 


1.571 


963 ± 144 


0.25±0.038 


0.034±0.005 



TABLE XV: Data for Case #3: equal-mass {q = 1) binary 
black holes with equal and opposite spins in the e*^^ — e'^' 
plane. Points A is from [l3 | while points -Bi — Bg are from 
[g|- Jtba/Lq is the ratio of the total radiated angular momen- 
tum to the initial orbital angular momentum, while E^^d/M 
is the ratio of the total radiated energy to the sum of the 
initial horizon masses. We do not include estimates of the 
radiated energy and angular momentum for point A because 
this simulation began at a different initial separation and or- 
bital angular momentum from points B\ — Bq. Herrmann et 
al. d claimed 15% errors on their reported numbers which 
we treat here as true Icr bars. 



APPENDIX B: SUPPLEMENTARY EQUATIONS 








|k|""'" 




-Brad /M 


7m ax 




0° 


854±128 


0.68±0.102 


0.046±0.0069 


192.3 


A2 


15° 


1401±210 


0.68±0.102 


0.044±0.0066 


189.5 


Ai 


30° 


2000±300 


0.67±0.101 


0.044±0.0066 


184.1 


Ai 


45° 


2030±305 


0.66±0.099 


0.043±0.0065 


177.3 


As 


60° 


1218±183 


0.65±0.098 


0.040±0.0060 


168.6 


Af, 


75° 


230±35 


0.64±0.096 


0.037±0.0056 


159.1 


A, 


90° 


1462±219 


0.62±0.093 


0.034±0.0051 


148.6 


As 


105° 


1979±297 


0.60±0.090 


0.033±0.0050 


138.6 


Aq 


120° 


1787±268 


0.58±0.087 


0.032±0.0048 


130.5 


Aw 


135° 


1234±185 


0.56±0.084 


0.030±0.0045 


124.1 


An 


150° 


689±103 


0.55±0.083 


0.029±0.0044 


119.5 


A12 


165° 


335 ±50 


0.55±0.083 


0.028±0.0042 


117.7 


Aw 


180° 


188±28 


0.55±0.083 


0.028±0.0042 


117.7 


Ai4 


195° 


157±24 


0.55±0.083 


0.028±0.0042 


120.5 


Ais 


210° 


173±26 


0.56±0.084 


0.030±0.0045 


125.5 


Aid 


225° 


223±33 


0.57±0.086 


0.032±0.0048 


132.7 


Air 


240° 


268±40 


0.59±0.089 


0.034±0.0051 


141.4 


Ais 


285° 


253±38 


0.65±0.098 


0.039±0.0059 


174.1 


Aig 


300° 


406±61 


0.66±0.099 


0.042±0.0063 


181.8 


A20 


315° 


399±60 


0.67±0.101 


0.045±0.0068 


187.7 


A21 


330° 


354±53 


0.68±0.102 


0.046±0.0069 


191.8 


A22 


345° 


459±69 


0.68±0.102 


0.046±0.0069 


193.2 



1. Relations between coefficients for Case #1 

Here are the relations between the "new" coefficients 
A, B, and C, and the original expansion coefficients 

1 mim2m3 \n1n2n3 

% : 



A-- 



0021000, 
^± I 

001|000| 



COS Q 



(Bla) 



TABLE XVL Data for Case #4: equal-mass {q = 1) BBHs 
belonging to the "S-Series" of 0|. The spins have magni- 
tudes a = 0.6 and orientations a = {—a, 0, 0) and b = 
(a sin cj>, 0, a cos 0). Jinai/M^ is the z— component of the final 
black hole's spin in units of the sum M of the initial hori- 
zon masses, and E^^g^lM is total radiated energy in units of 
M. Herrmann et al. [0] claimed 15% errors on their reported 
numbers which we treat here as true la bars. Tmax, measured 
in units of M, is an estimate of the merger time defined as 
the coordinate time between the beginning of the simulation 
and when the Newman-Penrose quantity ^'4 is maximized. 



B-- 



C-- 



001|000 , 003|000 
Ij^OOllOOOp 



1/lk' 



0021000, 



2 y |k°°^l°°"| 



sin2e(Blb) 



, 001|000 /, 002|001 , 003|000n 

^ ^ L + 2B 

I, 001 OOOi, ^ 



(Blc) 



, ^ . , , 1 u i 1 OOllOOO , , 0021000 

where B is the angle between k|_ ' and k|_ ' 







(i>a 


Ob 




Mf/M 






90° 
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0.94±0.0088 
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2. Relations between coefficients for Case #2 

Here are the relations between the "old" coefficients 
ymirnsmalnmsna ^nd the "new" Coefficients /(^'J) of sub- 



TABLE XVH: Data for Case #5: equal-mass {q = 
1) BBHs with generic spin orientations taken from 
[ll| . The initial spins have magnitudes a = 0.8 
and orientations given by traditional spherical coordi- 
nates, a = (a sin cos (/>a, a sin 6'a sin , a cos 6^^) and b = 
(a sin 6h cos 0(, , a sin 6^ sin (f>i, , a cos 9^, ) . 
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section IIVBI 



c j(0,0) _ 
cy(3,l), 



_ jOOO|000 

_ j200|000_|_ j020|000_ 1 ylOOllOO 
_ j200|000 _ y020|000_ 1_ ylOO|100_ 
_ jllO|000_ Jl00|010 

;2jioo|ooo 
;2joio|ooo 



1 j010|010 
1 ^010|010 



(B2a) 



i y300|000_ 3 j200|100_ 
i y020|100_ 1 Jll0|010 



.1): 



2 
I 
2 

3 J-030|000 
2J 2 



1 f210|000 
2^ 



iiolioo 



1 yl20l000 
1 j200|010 



3 j020|010 



(B2b) 



c j(3,3) _ 1 j300|000_ 1 j200|100_ 
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200 010 



3. Relations between coefficients for Case #3 

Here are the relations between the "old" coefficients 
f"nm.2 m,\mn 2n, ^^ic "ncw" Coefficients /(^'J) of sub- 
section |IVC] Cosine terms with even i,j appear in the 
expansion of scalars even under PX , the obscrvablcs m 
and S3. 



= y(0,0) _yOOO|000 

= y(2,0) _ j002|000^ j200|000_i j001|001_ 1 jlOO|100 {Qi) 
c j(2-2) _ J002|000_y200|000_ 1 jOOl |001 _|_ ]_ jlOO|100 



Cosine terms with odd i, j appear in the expansion of kj^ 
because it is a scalar odd under PX . 



4. Relations between coefficients for Case #4 



Here are the relations between the "old" coefficients 
jmima malnm ana and the "ucw" Coefficients /('■^' of sub- 
section llVDl Coefficients in the expansions of Jf^^,^\/M'^ 
and E-ci^^/M behave like those for m, which arc provided 
here. 
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We next provide expressions for the coefficients in the 
expansions of kj_, fc3, and |kp in Eq. in terms of the 
original coefficients of our general expansion (|25p . 
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Sine terms with odd i,j appear in the expansion of k^, 
because it is a pseudoscalar odd under PX. 
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because it is a pseudoscalar e?;e7i under PX . 
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APPENDIX C: THIRD-ORDER SPIN 
EXPANSIONS 



In Section [VBl wc identified 10 equal-mass initial spin 
configurations which when simulated could be used to 
calibrate all the coefiicients appearing in spin expansions 
of the 4 variables {«;, x, y, z} up to second order. Here we 
provide the corresponding third-order terms appearing in 
those same spin expansions. If desired, these formulae 
can be used to identify 12 additional equal-mass spin 
configurations with which these third-order terms may 
be calibrated. The third-order terms in the expansion 
for w {P = +1,X ^ +1) are 



-w'''\''\alb, + bla,)+w"''\''\alb, + bla,) 

yllllOOO/ , h h h \ , „„110|001 

-z«i°i|°i"(ai62«3 + &ia2&3)+«^'°"'°" 
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(ai&2^ + ^ia2a3) 
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yl01|100 
,002|001 



a-ibi{a 
0363 (a 



3^^3) 
3^^3) 



011|010 
003|000^„3 



-w—'—>a^b^{a:i + b^) 



(ai + bi). 



(Cla) 



The corresponding third-order terms in the expansion 
for a; (P = +1, X = —1) may be obtained from 
the above equation for w by making the substitution 

mim2m3|nin2n3 j, ^mim2m3\nin2n3 changing "-|-" 



W 



to " when it appears in parentheses: 
(...-...). 



•) 



The third-order terms in the expansion for y {P 
X = +1) are 
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,,020 
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''\a,bl + b,al) + y'''\''\a,bl + b,al) 

000(„3^^3)^y030|000(^3^^3)_ 



(Clb) 



The third order terms in the expansion for z (P = — 1, 
X ~ ~1) again may be obtained from the above equa- 
tion for y by making the substitution y™i™2m3|nin2n3 _^ 
^mim2Tn3|riin2"3 ^nd changing "+" to "-" when it ap- 
pears inside parentheses: (... + ...) ^ (... — .. .). 



APPENDIX D: MINIMUM- VARIANCE 
ESTIMATORS FOR THE SPIN COEFFICIENTS 

In SeclVlwe showed how a small number of simulations 
(10 in the equal-mass case, 16 for unequal masses) can 
be used to uniquely determine the 10 or 16 coefficients 
appearing to second order in the spin expansion. In the 
absence of systematic errors these are all the simulations 
that would be required, but further simulations may be 
useful once these errors are taken into account. In this 
Appendix, we will explicitly construct minimum- variance 
unbiased estimators for the coefficients in the spin expan- 
sion from N noisy simulations of known covariance. 

We proceed in two steps. In step one, we solve the 
problem under the assumption that the uncertainties in 
the initial spins are negligible compared to those in the fi- 
nal quantities. In step two. we explore how our approach 
might be modified to include the effects of these initial 
spin uncertainties. 



1. Neglecting initial spin uncertainties 

Imagine a generic final quantity / e {w, a;, y, z, u, w} 
such as those described in Section |Vl In the 
ith simulation, with initial spin configuration 

(») Ji) „(») M M M 



{al ,bi ,b2 ,b^ }, this quantity will have 

an estimated value /,; that when averaged over different 
assumptions for the systematic errors is equal to its true 
value ff. 



(Dl) 



The systematic errors introduce uncertainty in the values 
of / estimated for each simulation, and this uncertainty 
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can be described by the covariance matrix for /, 

F^J^{f^f,)~{f^){f,)■ (D2) 

Consider A'' initial spin configurations, which correspond 
to the A'' true final values fi [i = 1, . . . , N). If we trun- 
cate the spin expansion at finite order then, as seen in 
Eqs. ((55l) and ([58)) . these N values fi are a linear combi- 
nation of D spin expansion coefficients Cj . Thus we can 
write the relationship in matrix form as 

f = Ac (D3) 

where f is a column vector with N elements , c is a col- 
umn vector with D elements cj , and A is an x Z? ma- 
trix whose N rows consist of the D combinations of initial 
spin components in each simulation multiplying the co- 
efficients Cj. For example, when we rewrite Eq. (|55bp 
in the form of Eq. (jD3[) . wc have f = {xi, . . . , xn} and 
c = {xOOi|""o,2;"02|ooo^^^^^2.ioo|oiO|^ ^ith D = 6. The 

elements of the matrix A arc then easily read off from 

Eq. jUg, e.g. Ail = (4*^ " ^z^)- 

At least D simulations are required to determine the 
D coefHcients Cj . If this minimum number of simulations 
N = D are performed and the spin configurations of 
these simulations are chosen such that the D x D matrix 
A is invertible, then there is a unique estimator 

c = A^^f (D4) 

such that 

(c) = A-i(f) = A-^Ac = c. (D5) 

Here we have assumed that the initial spin components 
are known exactly and all the uncertainty lies in the es- 
timated final quantities f . We will relax this assumption 
later in this Appendix. The covariance matrix for 
this estimator is 

Cij — {ciCj) — {ci){cj) 

- Au!Aj,'{Mi)~AT^'Aj,\f,){f,) 

- ^iL'Aj'Fki , (D6) 

where here and throughout this Appendix wc adopt the 
Einstein convention of summing over repeated indices. 
In the absence of systematic errors (Fy = 0), D simula- 
tions would suffice to determine the spin coefficients with 
perfect accuracy {Cij = 0). 

With systematic errors present {Fij 7^ 0), a larger set 
of simulations {N > D) can be used to construct esti- 
mators c with lower variance provided these additional 
simulations arc at least partially uncorrelated. As the 
estimated final quantities f arc linear in the spin coeffi- 
cients, our estimator generalizes to 

c = Wf , (D7) 

where W is now a, D x N matrix of linear weights. For 
this estimator to be unbiased 



implying that W must satisfy the constraint 

WA = I , (D9) 

where I is the D x D identity matrix. The set of N spin 
configurations to be simulated must be chosen such that 
A has a left inverse. Eq. (|D9[) consists of constraints 
on the DN elements of W, leaving additional freedom 
for > to choose W to minimize the covariance 

= W,kW,i{hh) - W,kW,i{fk){fi) 

= W,kW,iFu, (DIO) 

or in matrix notation 

C = WFW'^. (Dll) 

As C is a real, symmetric matrix, it can be decomposed 
into a real, diagonal eigenvalue matrix cr and an orthog- 
onal eigenvector matrix O 

C = Oo-O^ . (D12) 

The columns of O give the uncorrelated linear combina- 
tions of estimators c^, while the elements of cr give the 
variances of these combinations. We specifically seek the 
weight matrix W that minimizes the sum of these eigen- 
values 

Trcr = TrC. (D13) 

We can determine this W by the method of Lagrange 
multipliers, with a Lagrangian given by 

L = Tr[C + A^(WA-I)], (D14) 

where A is a _D x _D matrix of Lagrange multipliers. Set- 
ting the partial derivatives dL/dXij to zero yields the 
D X D constraint equation (jD9[) . while dL/dWij = 
provides the additional D x N matrix equation 

2WF + AA^ = . (D15) 

Eqs. ([D9)) and (|DT5|) thus provide D{N + D) linear equa- 
tions for the elements of A and DN elements of W. 
Solving these equations, we find 

W = (A^F-iA)-iA^F-^ (D16) 

To summarize the analysis so far: If wc can neglect the 

(r) (r)-, 

errors in the initial spin components {al , 6, } that go 
into the construction of A, then the minimum variance 
unbiased estimators Ci for the spin-expansion coefficients 
Ci are given by Eq. (|D7[) . with weight matrix W given by 
Eq. (|D16[) . These optimal estimators c,; will have covari- 
ance matrix given by 



(c) = W(f) = WAc = c 



(D8) 



C = (A^F-iA)-^ (D17) 
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2. Including initial spin uncertainties 

Now let us consider the effect of errors in the initial 
spin components {a'f\b'f^}. One source of these errors 
is that numerical relativists do not know how to specify 
the proper initial data corresponding to physical binaries 
of arbitrary spin. This problem is usually addressed by 
waiting for the non-physical "junk radiation" to exit the 
system, after which the binary is presumed to settle down 
into a physical spin configuration. This physical spin con- 
figuration will generally be slightly different than the one 
relativists had intended to specify, introducing error into 
the initial spin components {ci[^\b[^^}. Techniques exist 
to measure BBH spins in simulations [1^ [53, [60, ^M, j 
so one could measure the initial spins after the junk ra- 
diation had exited the system, and use these components 
rather than those supposedly specified by the numerical 
initial data. However, while at large binary separations 
the initial spin components of the two black holes are well 
defined, at finite separations the components are gauge- 
dependent and different techniques for measuring them 
yield different results. 

To formally address the systematic errors in A we pro- 
mote it to an estimator A that on average will provide 
the correct values 

(iy) - Ay (D18) 

but will now have a non-zero covariance 

S,,ki = {A.jAki) - (iy)(ifci) . (D19) 

Deriving the truly optimal weight matrix W that si- 
multaneously minimizes contributions to the covariance 
from errors in both the final quantities f and initial spins 
A will be challenging. Since W is constructed from A as 
seen in Eq. (|D16[) . it is itself now an estimator W with 
its own covariance matrix 

T.jki = im^Wki) - (VKy) (Wki) . (D20) 

The covariance matrix of our estimator c will now be 
given by 



Cij — {ciCj) — {ci){cj) 

= (W^kWjJkfi) - c,Cj 



(D21) 



To make further progress, we make the possibly invalid 
assumption that errors in our estimators W and f are 
uncorrelated 

{mkW.ifkfi) = (mkWji) Cfkfi) ■ (D22) 

This allows us to reduce Eq. (jD21[) to 

= W.uWjiFki + T.kjihfi + T.kjiFki . (D23) 



The first term in Eq. (|D23p is the familiar error of 
Eq. (|D10[) , while the second and third terms proportional 
to Tijki reflect the increased covariance due to errors in 



the initial spins. One might next hope to insert this C 
into the Lagrangian L of Eq. (|D14p and obtain a new 
N X D matrix equation dL/dAij = to constrain A. 
Two problems immediately come to mind with this ap- 
proach. Firstly, once A has been promoted to an esti- 
mator Eqs. (|D9|) and (|D15|) become non-linear in W, A, 
and A and hence much more difficult to solve. Secondly, 
only 6A^ of the DN elements of A may be chosen in- 
dependently as there are only 6 initial spin components 
in each of the N simulations. Possibly this could be ad- 
dressed by adding new terms to the Lagrangian with new 
Lagrange multipliers, although this might be difficult to 
implement. 

Leaving the determination of a truly optimal estimator 
for future work, we choose to stick with the estimator c 
defined by the weight matrix W of Eq. (jD16p . This esti- 
mator should remain nearly optimal provided the initial 
spin errors are subdominant to the other errors coming 
from the simulations themselves and from truncating the 
spin expansion at finite order. Now we can account for 
the errors in the initial spins {af'\6^'"'} as follows. If 
these errors have a known probability distribution {e.g. 
if they are assumed to be Gaussian, with known covari- 
ance matrix), then we can Monte Carlo many realizations 
of {a^\l^^}, and use these to compute many realiza- 
tions of A. Next, by inserting these A-realizations into 
Eqs. (|D16[) and (|D7p . we obtain many realizations of c. 
The mean of these c-realizations is our best guess for c, 
while the covariance of these realizations gives an esti- 
mate of the "extra" uncertainty in the estimator c due 
to the initial spin uncertainties. We can add this "extra" 
covariance to the right hand side of Eq. (|D17p to estimate 
the total covariance Cij . 

As an alternative to Monte Carlo, we can make further 
analytic progress by assuming that the errors <5A in A 
are small, in which case the errors <5W in W can be 
linearized in SA. Defining 



N = A'^F^^A, 



we linearize Eq. (jD16[) to obtain 

= N-i(5A^ - ^A^F-^AN-^A^ 
-A^F-^AN^^A^IF^^ . 



(D24) 



(D25) 



This equation allows us_ to propagate errors and express 
the covariance Tij^i of W in terms of the covariance Sijki 
of A. Defining 



M = F-^AN-^A^ 



(D26) 



we obtain 



T^jki = {6W,,SWki) 

- 25cae/(F-lA),b(AN-l)d/ - 2ScaebMed 



■ 2SeafgM,,{F-'A) ft{AN-')dg] . (D27) 
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Inserting Eq. (|D27p in (jD23[) provides an approximate 
analytic expression for the "extra" covariance of the es- 
timator c caused by uncertainties in the initial spins. 

APPENDIX E: GENERALIZATION TO 
NON-CIRCULAR (ECCENTRIC) ORBITS 

We focused on circular orbits in the body of this paper 
as gravitational radiation is expected to circularize orbits 
of most astrophysical systems long before the final stage 
of the merger [43| . However, our approach readily gener- 
alizes to initially non-circular orbits so we felt that a few 
brief remarks on this subject might be appropriate here. 
Recall from Sec. that in the circular case, the initial 
conditions are specified by 8 dimcnsionlcss parameters: 
the mass ratio q and the dimcnsionlcss spins {a, b} at 
some initial instant labelled by ?/;• We can extend our 
spin expansion to non-circular orbits by specifying the 
difference in linear momentum p = — p^ between 
the two BBHs. As p lies in the orbital plane, our initial 
conditions arc now specified by 8-1-2=10 numbers. 

Apart from this modification, the analysis proceeds 
just as in Sec. |TT1 We define the same orthonormal triad 
{e'-^\e^^\e^'^^}, and consider the maps from the initial 
quantities to the final quantities 

f = f{^,p,,,q,a,,b,). (El) 



As in Sec. [Ill wc can constrain these maps through sym- 
metry considerations. Under parity P or exchange X, we 
have p ^ -p and {e^^\e^'^'>} -{e'l', e^^)}. Thus, the 
components Pi and P2 are invariant under both P and 
X. The maps therefore satisfy 

f{'>P,Pi,q,ai,b^) = (±)p/(^/>,p,;,g,a,,6,) (E2a) 

and 

f{^,Pi,q,ai,bi) = {±)pxf{ip,Pi, l/q,b^,ai). (E2b) 

Since the components Pi and P2 have eigenvalues 
of +1 under P and X, the series expansions intro- 
duced in Sec. IIIDI still hold, but now the coefficients 
ymimamahmans ^re functions of i/', q, and p^. It is proba- 
bly useful to Taylor expand these coefficinets around the 
point = where p^ ^.j^.^, is the linear momentum for 

a circular non-spinning orbit at orbital "separation" ijj. 
In the Newtonian limit, the three parameters {'0,Pi,P2} 
specify the semi-major axis, eccentricity, and longitude 
of pericenter associated with elliptical orbits. As the 
BBHs inspiral, they will trace a trajectory through this 
3-dimensional parameter space. Using this to relate coef- 
ficients defined at different points in the parameter space 
will be pursued in future work. 
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